|
| 1 | +## Term representation and simplification |
| 2 | + |
| 3 | +### Preliminary representation |
| 4 | + |
| 5 | +Expressions are stored as a tree where the nodes store their children in some efficient data structure. A node `n` is a non-leaf node `istree(n)` is true. For such nodes, `gethead` and `getargs` must be defined. These return the head and args of the term represented by the subtree `n`. |
| 6 | + |
| 7 | +A generic term is represented with the `Term` type. It simply holds the function `f` being called and a vector of arguments in the order they were passed to `f`. The operators `+` (and `-`), `*`, `/` and `^` create terms with special storage. |
| 8 | + |
| 9 | +Linear combinations such as $\alpha_1 x_1 + \alpha_2 x_2 +...+ \alpha_n x_n$ |
| 10 | +are represented by `Add(Dict(x_1 => \alpha_1, x_2 => \alpha_2, ..., x_n => \alpha_n))`. Now $x_n$ may themselves be other types mentioned in this section, but may not be an `Add`. When an `Add` is added to an `Add`, we merge their dictionaries and adding up matching coefficients. |
| 11 | + |
| 12 | +Similarly, $x_1^{m_1}x_2^{m_2}\ellipsisx_{m_n}$ is represented by |
| 13 | +`Mul(Dict(x_1 => m_1, x_2 => m_2,..., x_n => m_n))`. $x_i$ may not themselves be `Mul`, multiplying a Mul with another Mul returns a flattened Mul. |
| 14 | + |
| 15 | +$p / q$ is represented by `Div(p, q)`. The result of `*` on `Div` is maintainted as a `Div`. For example, `Div(p_1, q_1) * Div(p_2, q_2)` results in `Div(p_1 * p_2, q_1 * q_2)` and so on. The effect is, in `Div(p, q)`, `p` or `q` or, if they are Mul, any of their multiplicands is not a Div. So `Mul`s must always be nested inside a `Div` and can never show up immediately wrapping it. |
| 16 | + |
| 17 | + |
| 18 | +### Polynomial representation |
| 19 | + |
| 20 | +Packages like DynamicPolynomials.jl provide even more efficient representations for polynomials. They also have efficient operations such as multi-variate polynomial GCD. However, DynamicPolynomials can only represent flat polynomials. For example, `(x-3)*(x+5)` can only be represented as `(x^2) + 15 - 8x`. Moreover, DynamicPolynomials does not have ways to represent generic Terms such as `sin(x-y)` in the tree. |
| 21 | + |
| 22 | +To reconcile these differences while being able to use the efficient representation of DynamicPolynomials we have the `PolyForm` type. This type holds a polynomial and the mappings necessary to present them as a SymbolicUtils expression. The mappings constructed for the conversion are 1) a bijection from DynamicPolynomials PolyVar type to a Symbolics `Sym`, and 2) a mapping from `Sym`s to non-polynomial terms that the `Sym`s stand-in for. These terms may themselves contain PolyForm if there are polynomials inside them. The mappings are transiently global, that is, when all references to the mappings go out of scope, they are released and re-created. |
| 23 | + |
| 24 | +```julia |
| 25 | +julia> @syms x y |
| 26 | +(x, y) |
| 27 | + |
| 28 | +julia> PolyForm((x-3)*(y-5)) |
| 29 | +x*y + 15 - 5x - 3y |
| 30 | +``` |
| 31 | + |
| 32 | +Further, generic terms such as `sin` or `cos` terms are replaced with a symbol when representing the polynomial. |
| 33 | + |
| 34 | +```julia |
| 35 | +julia> p = PolyForm((sin(x) + cos(x))^2) |
| 36 | +(cos(x)^2) + 2cos(x)*sin(x) + (sin(x)^2) |
| 37 | + |
| 38 | +julia> p.p # this is the actual DynamicPolynomial stored |
| 39 | +cos_3658410937268741549² + 2cos_3658410937268741549sin_10720964503106793468 + sin_10720964503106793468² |
| 40 | +``` |
| 41 | + |
| 42 | +By default, polynomials inside generic terms are not also converted to PolyForm. For example, |
| 43 | + |
| 44 | +```julia |
| 45 | +julia> PolyForm(sin((x-3)*(y-5))) |
| 46 | +sin((x - 3)*(y - 5)) |
| 47 | +``` |
| 48 | +But you can pass in the `recurse=true` keyword argument to get this going. |
| 49 | + |
| 50 | +```julia |
| 51 | +julia> PolyForm(sin((x-3)*(y-5)), recurse=true) |
| 52 | +sin(x*y + 15 - 5x - 3y) |
| 53 | +``` |
| 54 | + |
| 55 | +Polynomials are constructed by first turning symbols and generic terms into Polynomial variables and then applying the `+`, `*`, `^` operations on these variables. You can control which of these operations are actually constructed with the `Fs` keyword argument. It is a Union type of the functions to apply. For example, let's say you want to turn terms into polynomials by only applying the `+` and `^` operations, and want to preserve `*` operations as-is, you could pass in `Fs=Union{typeof(+), typeof(^)}` |
| 56 | + |
| 57 | +```julia |
| 58 | +julia> PolyForm((x+y)^2*(x-y), Fs=Union{typeof(+), typeof(^)}, recurse=true) |
| 59 | +((x - (y))*((x^2) + 2x*y + (y^2))) |
| 60 | +``` |
| 61 | + |
| 62 | +Note that in this case `recurse=true` was necessary as otherwise the polynomialization would not descend into the `*` operation as it is now considered a generic operation. |
| 63 | + |
| 64 | +### Simplifying fractions |
| 65 | + |
| 66 | +`simplify_fractions(expr)` recurses through `expr` finding `Div`s and simplifying them using polynomial divison. |
| 67 | + |
| 68 | +First the factors of the numerators and the denominators are converted into PolyForm objects, then numerators and denominators are divided by their respective pairwise GCDs. The conversion of the numerator and denominator into PolyForm is set up so that `simplify_fractions` does not result in increase in the expression size due to polynomial expansion. Specifically, the factors are individually converted into PolyForm objects, and any powers of polynomial is not expanded, but the divison process repeatedly divides them as many times as the power. |
| 69 | + |
| 70 | + |
| 71 | +```julia |
| 72 | +julia> simplify_fractions((x*y+5x+3y+15)/((x+3)*(x-4))) |
| 73 | +(5.0 + y) / (x - 4) |
| 74 | + |
| 75 | +julia> simplify_fractions((x*y+5x+3y+15)^2/((x+3)*(x-4))) |
| 76 | +((5.0 + y)*(15 + 5x + x*y + 3y)) / (x - 4) |
| 77 | + |
| 78 | +julia> simplify_fractions(3/(x+3) + x/(x+3)) |
| 79 | +1 |
| 80 | + |
| 81 | +julia> simplify_fractions(sin(x)/cos(x) + cos(x)/sin(x)) |
| 82 | +(cos(x)^2 + sin(x)^2) / (cos(x)*sin(x)) |
| 83 | + |
| 84 | +julia> simplify(ans) |
| 85 | +1 / (cos(x)*sin(x)) |
| 86 | +``` |
0 commit comments