Movimentum - Polynomial folding 3

5. Plus

Using "ASCII art," here is a table that shows how addition must work. Again, the goal is to merge polynomial subexpressions of a single variable at the left, and collect everything else in a single "rest" expression on the right. Be aware that E and F might also be polynomials, if their variables are different from V!

    //        | Q[V]+F              Q[V]            F
    // -------+------------------------------------------------
    // P[V]+E | {P[V]+Q[V]}+(E+F)   {P[V]+Q[V]}+E   P[V]+(E+F)
    //        |
    // P[V]   | {P[V]+Q[V]}+F       {P[V]+Q[V]}     P[V]+E
    //        |
    // E      | Q[V]+(E+F)          Q[V]+E          E+F

Why don't I consider entries with unary minus, or even square roots?

• Square roots, of course, cannot be merged together by addition, so handling them separately will not create smaller, more compact expressions.
• For unary minus, the argument is a little trickier: First, there cannot be a –P[V] subexpression anywhere—it would have been rewritten to {–P[V]}, i.e., a simple polynomial with inverted coefficients, on a lower level. The E's and F's could be expressions with a leading unary minus, but moving their minus upwards will not reduce the complexity of the expression in general. In specific cases, such a rewrite might get rid of expressions—e.g., if E happens to be –F or the like. But we do not optimize such "rare occurrences" here.

Polynomial folding is a heuristics. Like constant folding, it tries to standardize expressions quickly. However, if it oversees certain simplifications, this is not a big problem, as long as we do not rely on simplification guarantees.

Coding the table for the plus operator is easy—for example, here are the three entries for the first line of the table:

    new StandardExpressionRewrite("(P+E)+(Q+F)", _plusRewrites,
        (p + e) + (q + f),
        HaveSameVar(p, q),
        m => ((m & p) + (m & q)) + ((m & e) + (m & f)));
    //                PO         AE         AE

    new StandardExpressionRewrite("(P+E)+Q", _plusRewrites,
        (p + e) + q,
        HaveSameVar(p, q),
        m => ((m & p) + (m & q)) + (m & e));

    new StandardExpressionRewrite("(P+E)+F", _plusRewrites,
        (p + e) + f,
        m => (m & p) + ((m & e) + (m & f)));
    }

Below the first rewrite rule, I have marked the + operators with their "type":

• PO is a polynomial operator, i.e., it will merge its two operands into a single polynomial.
• AE, on the other hand, is an AbstractExpr operator, i.e., it will create an expression tree (a BinaryExpression).

The actual visiting code is similar to the square and unary minus operators:

    public IAbstractExpr Visit(Plus op, IAbstractExpr lhs,
                               IAbstractExpr rhs, Ignore p) {
        return Rewrite(new BinaryExpression(lhs, op, rhs),
                       "Cannot handle " + op, _plusRewrites);
    }

6. Times

My first attempt for multiplication was the same as for addition: Create a table that contains all combinations of "interesting expression structures," and explicitly define the result. From the square operator, we learn that we should handle square roots separately if they occur on both sides of the operator. In that first attempt, I identified 10 subexpression forms that I needed to distinguish: P+√R, P+–√R, P+√E, P+–√E, P+E, P+–E, P, C, –E, and general E. This gives us a table of 100 entries (or even more, because we have to distinguish (P+√R)·(Q+√R) from (P+√R)·(Q+√S), where all of P, Q, R, S are polynomials of the same variable). This is too much work—specifying, coding, and last, but not least, testing.

I thought a little bit and came up with the idea of stage-wise processing—it is what we do manually to "simplify" such expressions. For example, if we have (x + √4x) · (x - √x), we will

• first multiply the two terms which gives us a sum of four terms;
• then simplify each term;
• and finally the whole sum.

The result from the first step can be a sum of up to four terms, i.e., an expression tree of depth at most 3. Various factors in the sums might need simplification—e.g. in our example, where √4x · √x will have to be merged to √4x², which means that we have to run the tree through folding again at least three levels down from the current root. Instead of introducing a special mechanism for this, we run the visitor over the complete result—this will uselessly visit expressions below depth 3, but we will have to optimize away useless visits anyway somewhen in the future, so we will ignore this overhead right now.

The three steps above are only a rough plan. Actually, we need to consider a little bit more what "simplifications" we have to apply to the result of the first step. Here is my suggestion for a three-step algorithm.

In the first version of this posting, I bragged "similarly to all previous algorithms, I do not give any proof why this accomplishes what we want; it is "more or less obvious," ...". Well, in addition to being obvious, it was also wrong. The tables for Step III claimed, among others, that P√Q could be transformed to √(P²Q). However, this cannot be true: The latter is never negative (as we define sqrt to be the positive square root), but the former can of course be negative—for example, if the polynomial P is x – 1, and x happens to be zero. Here is a better table; but still no proof—I risk keeping this a mere software project, not a scientific undertaking, although I'd like to do all these proofs. Maybe I find time for them somewhen later ...

    // Step I: Distribute.
    //
    //          Q+F                 F
    //
    //  P+E     {PQ}+Q*E+P*F+E*F    P*F+E*F
    //
    //  E       Q*E+E*F             =
    //

    // Step II: Lift unary minus to top.
    //
    //          -F       F
    //
    //  -E      E*F      -(E*F)
    //
    //  E       -(E*F)   =
    //

    // Step III: Pull up square roots as far as possible.
    //
    //          √Q              √F              Q       D               F
    //
    //  √P      P=Q:  P         √(P*F)          =       D<0: -√{PD²}    =
    //          else: √{PQ}                             else: √{PD²}
    //
    //  √E      √(Q*E)          E=F:  E         =       D<0: -√({D²}*E) =
    //                          else: √(E*F)            else: √({D²}*E)
    //
    //  P       =               =               {P*Q}   {P*D}           =
    //                                                X
    //  C       C<0: -√({C²Q}   C<0: -√({C²}*F) {C*Q}   {C*D}           =  
    //          else: √({C²Q}   else: √({C²}*F)
    //
    //  E       =             =                 =       =               =
    //

An equals sign, in the tables above, simple means that no rewrite is done for the indicated product. Together, these three steps require 3 + 3 + 15 = 21 rules—a lot less than the 100 or more needed in my first attempt!
The four rewritings around the X can be handled by the same rule, which simply matches the product of two polynomials—also constants are polynomials!

Here is the code for Step I: The rules are almost like earlier ones—but in addition, there is a call to RecursivelyFold on the result:

    new StandardExpressionRewrite("(P+E)*(Q+F)", _timesRewrites,
        (p + e) * (q + f),
        HaveSameVar(p, q),
        m => RecursivelyFold((m & p) * (m & q)
                            + (m & q) * (m & e)
                            + (m & p) * (m & f)
                            + (m & e) * (m & f)));
    new StandardExpressionRewrite("(P+E)*F", _timesRewrites,
        (p + e) * f,
        m => RecursivelyFold((m & p) * (m & f)
                            + (m & e) * (m & f)));
    new StandardExpressionRewrite("E*(Q+F)", _timesRewrites,
        e * (q + f),
        m => RecursivelyFold((m & q) * (m & e)
                            + (m & e) * (m & f)));

RecursivelyFold is a single line method:

    private static readonly PolynomialFoldingVisitor
        _recursiveFolder = new PolynomialFoldingVisitor();

    private static IAbstractExpr RecursivelyFold(AbstractExpr expr) {
        return expr.Accept(_recursiveFolder);
    }

Here is an excerpt from Step III:

    new StandardExpressionRewrite("√P*√P", _timesRewrites,
        sqrtP * sqrtP,
        m => RecursivelyFold(m & p));
    new StandardExpressionRewrite("√P*√Q", _timesRewrites,
        sqrtP * sqrtQ,
        HaveSameVar(p, q),
        m => PosSqrtAndRecursivelyFold((m & p) * (m & q)));
    new StandardExpressionRewrite("√P*√F", _timesRewrites,
        sqrtP * sqrtF,
        m => PosSqrtAndRecursivelyFold((m & p) * (m & f)));
    // √P*Q --> =
    new StandardExpressionRewrite("√P*D,D<0", _timesRewrites,
        sqrtP * d, m => (m & d).Value < 0,
        m => -PosSqrtAndRecursivelyFold(
                    (m & d).Value * (m & d).Value * (m & p)
        ));

PosSqrtAndRecursivelyFold wraps the parameter in a PositiveSquareroot after recursively folding it:

    private static IAbstractExpr PosSqrtAndRecursivelyFold(AbstractExpr expr) {
        return new UnaryExpression(RecursivelyFold(expr), new PositiveSquareroot());
    }

If you look closely, you will notice that I put all the rules for multiplication into a single list, called _timesRewrites. This means that later stages, i.e. the RecursiveFold calls, will run over the rules for earlier stages again. This is actually useless and could be optimized away, but I guess that that would require even more complex algorithms than the ones I invent here. So I leave it as it is.

Of course, we also need an implementation of the Visit method: But it is similar to the one for the + operator.

Together with more than one hundred unit tests, this code rewrites binary expressions with a Times operator to our intended P+E format!

6. Divide

Division, finally, is very simple:

• If the division is by a constant, we replace it with a multiplication of the reciprocal value.
• Otherwise, we leave the expression alone (we do not try polynomial division or the like!):

    new StandardExpressionRewrite("E/C", _divideRewrites,
        e / c,
        m => RecursivelyFold(
              (m & e) * Polynomial.CreateConstant(1 / (m & c).Value)
            )
    );
    new StandardExpressionRewrite("E", _divideRewrites,
        e,
        m => m & e);

The corresponding Visit method is, of course, like the ones for plus and times.

And this concludes the visitor for rewriting polynomials.