Movimentum - Introducing polynomials

In the previous posting, we ended up with six unsolvable constraints that can be condensed to the following:

! {EqualsZeroConstraint} 0 = f₁(x)
! {EqualsZeroConstraint} 0 = f₂(x)
! {EqualsZeroConstraint} 0 = f₃(x)
! {AtLeastZeroConstraint} 0 <= f₄(x)
! {AtLeastZeroConstraint} 0 <= f₅(x)
! {AtLeastZeroConstraint} 0 <= f₆(x)

The f_i are expressions that contain only the variable x and a few operators and constants—nothing else. Now, at least for the first three constraints, this looks definitely like a textbook problem: We have all learned various "iteration" or "numerical" algorithms to find zeros of a function—at least "Newton", bisection, secant come to mind. And therefore I tried it ... which has cost me a few days of my life, because all these methods do not work (at least when implemented in a straightforward manner). The reason can be seen from a plot of one of the functions, which I made with wolframalpha.com:

The presence of square roots leads to functions that can have

• vertical tangents at some points;
• are defined only for limited intervals of x;

Moreover, many functions are even—this is bad for algorithms that need two function values of opposite sign to start their search, like bisection. You can try all your skills to find a general numerical solving method that finds the zeros for such functions—I gave up.

(By the way, it is interesting that wolframalpha can find out that the function is defined on the real numbers only if x is at least 100; but still, it gives 99.99999999986 as a solution. I think this means that they do use a numerical algorithm!).

Instead of finding a numerical algorithm, I sticked to the algebraic approach, i.e., formal rewriting of constraints. However, with our current constraints, this is virtually impossible. The reason is that the same constraint can be formulated in so many different ways that we cannot hope to write a manageable, sufficiently small number of rules to handle all important cases. Let us look at this problem in concrete at the first unsolvable constraint from the last posting:

0=200+-(--(200+10*x*0.923879532511287+0)+sqrt(10000+-((200+---(200+10*x*0.38268343236509+0))²+0)))

At first, this is a chaotic convolution of parentheses, numbers, and a few xs. But of course, with some simple, if a little tedious, simplifications this can be rewritten as

0 = 9.23879532511287 x + sqrt(10000 – 3.8268343236509² x²)

or, in abstract form,

0 = P[x] + sqrt(Q[x])

where P and Q are two terms containing only the variable x.

Now, a constraint of this structure can easily be rewritten to a squareroot-free form (there is a little bit missing here which we will handle later):

–P[x] = sqrt(Q[x])
P[x]² = Q[x]
0 = Q[x] – P[x]²

And from here, we—i.e., our solver—can explore new ways of solving, e.g. solving quadratic equations (if P and Q are simple enough).

But—big "but": In order to see a general pattern like '0 = P[x] + sqrt(Q[x])' in a constraint, we must standardize expressions much more than we did up to now.

But (again) we are fortunate—we do not have to look too far for one sort of standardizable expression that may help us: Polynomials. A polynomial can combine many of our current operators in a single, normalized way.

("Normalized" is even better than "standardized:" The latter means "can mostly be brought to one of a a handful of similar form." The former means "can always be brought to a single form.")

However, there are two problems if we want to standardize (or even normalize) our expressions to polynomials:

• First, we have expressions that cannot be converted to polynomials—e.g., square roots!
• Second, even rewriting to polynomials is complicated if there is more than one variable involved.

Regarding the first problem, my idea is to write each expression to the following:

E → P + E₁ + E₂ + ...

where P is a polynomial, and the E_i are the "other expressions." Recursively, I want each subexpression of E to have the same standard form.

Regarding the second item, we simply limit our scope to polynomials with single variables. All other expressions are, at least for the moment, handled as, well, "other expressions."

Polynomials

What do we need?

• A class that represents a polynomial.
• A way to include variables and constants into the design—after all, a single variable and a constant are (quite degenerate) polynomials, too.
• A simplification procedure that reduces general expression to the P + E_i form.

The following code shows our simple polynomial class for polynomials of a single variable. The _variable points to the variable. _coefficents contains the coefficients in "mathematical order," i.e., _coefficents[0] contains the highest-order coefficient, _coefficents[1] the next-lower one, etc. For example, the polynomial x³ + 2x² – 3x + 100 stores it coefficients in an array with

• _coefficients[0] = 1.0
• _coefficients[1] = 2.0
• _coefficients[2] = –3.0
• _coefficients[3] = 100.0

Likewise, x² – 3.14 is stored as [1.0, 0.0, -3.14]. Finally, the degree is stored in a separate field because I need it at two or three places and did not want to write _coefficients.Length - 1 more than once:

private class Polynomial : AbstractExpr {
    private readonly Variable _variable;
    private readonly double[] _coefficients;
    private readonly int _degree;
    internal GeneralPolynomial(Variable variable,
            IEnumerable<double> coefficients) {
        _coefficients = coefficients.ToArray();
        if (_coefficients[0].Near(0)) {
            throw new ArgumentException("Top coefficient must not be 0");
        }
        _degree = _coefficients.Length - 1;
        _variable = variable;
    }

    public override Variable Var { get { return _variable; } }

    public override int Degree { get { return _degree; } }

    public override double Coefficient(int power) {
        return _coefficients[_degree - power];
    }

    public override bool Equals(object obj) { ... }

    public override int GetHashCode() { ... }

    public override IEnumerable<double> Coefficients {
        get { return _coefficients; }
    }

    public override TResult Accept<TParameter, TResult>(
            ISolverModelExprVisitor<TParameter, TResult> visitor,
            TParameter p) {
        return visitor.Visit(this, p);
    }
}

The Accept method creates a new problem: We have to add a new method to the visitor interface ...

    TResult Visit(Polynomial polynomial, TParameter parameter);

... and consequently, we have to update all our existing visitors. This is a well-known effect of the visitor pattern: Modifications of the visited class structure is disruptive. However, in closed scenarios like ours, this is actually a benefit: It forces us to think about all existing operations for the new Polynomial class. Here are all the visitors that we have to update:

• ConstantFoldingVisitor: If we happen to encounter a polynomial of degree zero, we convert it to a constant. Otherwise, it can never be folded to a constant, because it contains a variable:

    public AbstractExpr Visit(Polynomial polynomial, Ignore p) {
        return polynomial.Degree == 0
            ? new Constant(polynomial.Coefficient(0))
            : polynomial;
    }

• EvaluationVisitor: The polynomial is evaluated using Horner's rule:

    public double Visit(Polynomial polynomial, Ignore parameter) {
        // Evaluate by Horner's rule.
        double value = GetValue(polynomial.Var);
        double result = polynomial.Coefficient(polynomial.Degree);
        for (int i = polynomial.Degree - 1; i >= 0; i--) {
            result = value * result + polynomial.Coefficient(i);
        }
        return result;
    }

• FindFormalSquarerootVisitor: A polynomial never contains a formal square root:

    public FindFormalSquarerootVisitor Visit(Polynomial poly, Ignore p) {
        return this;
    }

• RewritingVisitor: The rewriting result is exactly the same as for evaluation, only that the * and + operators do not create double numbers, but expressions. (The following code is inefficient when the variable is rewritten as a constant: In that case, a complete expression tree is created, which will later be folded to a constant by the ConstantFoldingVisitor. This case should, therefore, be handled separately by an EvaluationVisitor):

    public AbstractExpr Visit(Polynomial polynomial, Ignore p) {
        if (polynomial.Var.Equals(_from)) {
            // Evaluate by Horner's rule.
            AbstractExpr result =
                new Constant(polynomial.Coefficient(polynomial.Degree));
            for (int i = polynomial.Degree - 1; i >= 0; i--) {
                result = _to * result
                       + new Constant(polynomial.Coefficient(i));
            }
            return result;
        } else {
            return polynomial;
        }
    }

• ToStringVisitor: Converting a polynomial to a more or less readable form takes about 50 lines of code that I do not show here. My goal was threefold:
1. Create readable output—therefore, there is e.g. no x^1 in the output, but a plain x; and terms with zero coefficients are suppressed;
2. Create output that can directly be pasted into wolframalpha.com for evaluation and solving—this has helped me a number of times to find out whether some constraint, after rewriting, still had the desired solution;
3. Indicate that there is a Polynomial object in an expression tree. For this, I use (( and )) as enclosing parentheses. I could not use brackets or curly braces because of the previous item, but I hope that the doubled parentheses will stand out enough.
   
• VariableDegreeVisitor: Finally, this visitor simply returns a single-valued dictionary with the polynomial's variable and a suitable degree—not really interesting.

This was easy. The next step is not: Integrating the new class with existing classes that already are a sort of polynomials: Constants and variables.