# Polynomial evaluation techniques that have been developed

algebra c# expression-trees linq math

### Question

I'm trying to come up with an elegant way to handle some generated polynomials. Here's the situation we'll focus on (exclusively) for this question:

1. order is a parameter in generating an nth order polynomial, where n:=order + 1.
2. i is an integer parameter in the range 0..n
3. The polynomial has zeros at x_j, where j = 1..n and j â‰  i (it should be clear at this point that StackOverflow needs a new feature or it's present and I don't know it)
4. The polynomial evaluates to 1 at x_i.

Since this particular code example generates x_1 .. x_n, I'll explain how they're found in the code. The points are evenly spaced `x_j = j * elementSize / order` apart, where `n = order + 1`.

I generate a `Func<double, double>` to evaluate this polynomialÂ¹.

``````private static Func<double, double> GeneratePsi(double elementSize, int order, int i)
{
if (order < 1)
throw new ArgumentOutOfRangeException("order", "order must be greater than 0.");

if (i < 0)
throw new ArgumentOutOfRangeException("i", "i cannot be less than zero.");
if (i > order)
throw new ArgumentException("i", "i cannot be greater than order");

ParameterExpression xp = Expression.Parameter(typeof(double), "x");

// generate the terms of the factored polynomial in form (x_j - x)
List<Expression> factors = new List<Expression>();
for (int j = 0; j <= order; j++)
{
if (j == i)
continue;

double p = j * elementSize / order;
}

// evaluate the result at the point x_i to get scaleInv=1.0/scale.
double xi = i * elementSize / order;
double scaleInv = Enumerable.Range(0, order + 1).Aggregate(0.0, (product, j) => product * (j == i ? 1.0 : (j * elementSize / order - xi)));

/* generate an expression to evaluate
*   (x_0 - x) * (x_1 - x) .. (x_n - x) / (x_i - x)
* obviously the term (x_i - x) is cancelled in this result, but included here to make the result clear
*/
Expression expr = factors.Skip(1).Aggregate(factors[0], Expression.Multiply);
// multiplying by scale forces the condition f(x_i)=1
expr = Expression.Multiply(Expression.Constant(1.0 / scaleInv), expr);

Expression<Func<double, double>> lambdaMethod = Expression.Lambda<Func<double, double>>(expr, xp);
return lambdaMethod.Compile();
}
``````

The problem: I also need to evaluate Ïˆâ€²=dÏˆ/dx. To do this, I can rewrite Ïˆ=scaleÃ—(x_0 - x)(x_1 - x)Ã—..Ã—(x_n - x)/(x_i - x) in the form Ïˆ=Î±_nÃ—x^n + Î±_nÃ—x^(n-1) + .. + Î±_1Ã—x + Î±_0. This gives Ïˆâ€²=nÃ—Î±_nÃ—x^(n-1) + (n-1)Ã—Î±_nÃ—x^(n-2) + .. + 1Ã—Î±_1.

For computational reasons, we can rewrite the final answer without calls to `Math.Pow` by writing Ïˆâ€²=xÃ—(xÃ—(xÃ—(..) - Î²_2) - Î²_1) - Î²_0.

To do all this "trickery" (all very basic algebra), I need a clean way to:

1. Expand a factored `Expression` containing `ConstantExpression` and `ParameterExpression` leaves and basic mathematical operations (end up either `BinaryExpression` with the `NodeType` set to the operation) - the result here can include `InvocationExpression` elements to the `MethodInfo` for `Math.Pow` which we'll handle in a special manner throughout.
2. Then I take the derivative with respect to some specified `ParameterExpression`. Terms in the result where the right hand side parameter to an invocation of `Math.Pow` was the constant 2 are replaced by the `ConstantExpression(2)` multiplied by what was the left hand side (the invocation of `Math.Pow(x,1)` is removed). Terms in the result that become zero because they were constant with respect to x are removed.
3. Then factor out the instances of some specific `ParameterExpression` where they occur as the left hand side parameter to an invocation of `Math.Pow`. When the right hand side of the invocation becomes a `ConstantExpression` with the value `1`, we replace the invocation with just the `ParameterExpression` itself.

Â¹ In the future, I'd like the method to take a `ParameterExpression` and return an `Expression` that evaluates based on that parameter. That way I can aggregate generated functions. I'm not there yet. Â² In the future, I hope to release a general library for working with LINQ Expressions as symbolic math.

1
5
10/27/2009 2:54:03 PM

#### Fastest Entity Framework Extensions

I wrote the basics of several symbolic math features using the ExpressionVisitor type in .NET 4. It's not perfect, but it looks like the foundation of a viable solution.

• `Symbolic` is a public static class exposing methods like `Expand`, `Simplify`, and `PartialDerivative`
• `ExpandVisitor` is an internal helper type that expands expressions
• `SimplifyVisitor` is an internal helper type that simplifies expressions
• `DerivativeVisitor` is an internal helper type that takes the derivative of an expression
• `ListPrintVisitor` is an internal helper type that converts an `Expression` to a prefix notation with a Lisp syntax

## `Symbolic`

``````public static class Symbolic
{
public static Expression Expand(Expression expression)
{
return new ExpandVisitor().Visit(expression);
}

public static Expression Simplify(Expression expression)
{
return new SimplifyVisitor().Visit(expression);
}

public static Expression PartialDerivative(Expression expression, ParameterExpression parameter)
{
bool totalDerivative = false;
return new DerivativeVisitor(parameter, totalDerivative).Visit(expression);
}

public static string ToString(Expression expression)
{
ConstantExpression result = (ConstantExpression)new ListPrintVisitor().Visit(expression);
return result.Value.ToString();
}
}
``````

## Expanding expressions with `ExpandVisitor`

``````internal class ExpandVisitor : ExpressionVisitor
{
protected override Expression VisitBinary(BinaryExpression node)
{
var left = Visit(node.Left);
var right = Visit(node.Right);

if (node.NodeType == ExpressionType.Multiply)
{
Expression[] leftNodes = GetAddedNodes(left).ToArray();
Expression[] rightNodes = GetAddedNodes(right).ToArray();
var result =
leftNodes
.SelectMany(x => rightNodes.Select(y => Expression.Multiply(x, y)))
.Aggregate((sum, term) => Expression.Add(sum, term));

return result;
}

if (node.Left == left && node.Right == right)
return node;

return Expression.MakeBinary(node.NodeType, left, right, node.IsLiftedToNull, node.Method, node.Conversion);
}

/// <summary>
/// Treats the <paramref name="node"/> as the sum (or difference) of one or more child nodes and returns the
/// the individual addends in the sum.
/// </summary>
private static IEnumerable<Expression> GetAddedNodes(Expression node)
{
BinaryExpression binary = node as BinaryExpression;
if (binary != null)
{
switch (binary.NodeType)
{
foreach (var n in GetAddedNodes(binary.Left))
yield return n;

foreach (var n in GetAddedNodes(binary.Right))
yield return n;

yield break;

case ExpressionType.Subtract:
foreach (var n in GetAddedNodes(binary.Left))
yield return n;

foreach (var n in GetAddedNodes(binary.Right))
yield return Expression.Negate(n);

yield break;

default:
break;
}
}

yield return node;
}
}
``````

## Taking a derivative with `DerivativeVisitor`

``````internal class DerivativeVisitor : ExpressionVisitor
{
private ParameterExpression _parameter;
private bool _totalDerivative;

public DerivativeVisitor(ParameterExpression parameter, bool totalDerivative)
{
if (_totalDerivative)
throw new NotImplementedException();

_parameter = parameter;
_totalDerivative = totalDerivative;
}

protected override Expression VisitBinary(BinaryExpression node)
{
switch (node.NodeType)
{
case ExpressionType.Subtract:
return Expression.MakeBinary(node.NodeType, Visit(node.Left), Visit(node.Right));

case ExpressionType.Multiply:
return Expression.Add(Expression.Multiply(node.Left, Visit(node.Right)), Expression.Multiply(Visit(node.Left), node.Right));

case ExpressionType.Divide:
return Expression.Divide(Expression.Subtract(Expression.Multiply(Visit(node.Left), node.Right), Expression.Multiply(node.Left, Visit(node.Right))), Expression.Power(node.Right, Expression.Constant(2)));

case ExpressionType.Power:
if (node.Right is ConstantExpression)
{
return Expression.Multiply(node.Right, Expression.Multiply(Visit(node.Left), Expression.Subtract(node.Right, Expression.Constant(1))));
}
else if (node.Left is ConstantExpression)
{
return Expression.Multiply(node, MathExpressions.Log(node.Left));
}
else
{
Expression.Multiply(Visit(node.Left), Expression.Divide(node.Right, node.Left)),
Expression.Multiply(Visit(node.Right), MathExpressions.Log(node.Left))
));
}

default:
throw new NotImplementedException();
}
}

protected override Expression VisitConstant(ConstantExpression node)
{
return MathExpressions.Zero;
}

protected override Expression VisitInvocation(InvocationExpression node)
{
MemberExpression memberExpression = node.Expression as MemberExpression;
if (memberExpression != null)
{
var member = memberExpression.Member;
if (member.DeclaringType != typeof(Math))
throw new NotImplementedException();

switch (member.Name)
{
case "Log":
return Expression.Divide(Visit(node.Expression), node.Expression);

case "Log10":
return Expression.Divide(Visit(node.Expression), Expression.Multiply(Expression.Constant(Math.Log(10)), node.Expression));

case "Exp":
case "Sin":
case "Cos":
default:
throw new NotImplementedException();
}
}

throw new NotImplementedException();
}

protected override Expression VisitParameter(ParameterExpression node)
{
if (node == _parameter)
return MathExpressions.One;

return MathExpressions.Zero;
}
}
``````

## Simplifying expressions with `SimplifyVisitor`

``````internal class SimplifyVisitor : ExpressionVisitor
{
protected override Expression VisitBinary(BinaryExpression node)
{
var left = Visit(node.Left);
var right = Visit(node.Right);

ConstantExpression leftConstant = left as ConstantExpression;
ConstantExpression rightConstant = right as ConstantExpression;
if (leftConstant != null && rightConstant != null
&& (leftConstant.Value is double) && (rightConstant.Value is double))
{
double leftValue = (double)leftConstant.Value;
double rightValue = (double)rightConstant.Value;

switch (node.NodeType)
{
return Expression.Constant(leftValue + rightValue);
case ExpressionType.Subtract:
return Expression.Constant(leftValue - rightValue);
case ExpressionType.Multiply:
return Expression.Constant(leftValue * rightValue);
case ExpressionType.Divide:
return Expression.Constant(leftValue / rightValue);
default:
throw new NotImplementedException();
}
}

switch (node.NodeType)
{
if (IsZero(left))
return right;
if (IsZero(right))
return left;
break;

case ExpressionType.Subtract:
if (IsZero(left))
return Expression.Negate(right);
if (IsZero(right))
return left;
break;

case ExpressionType.Multiply:
if (IsZero(left) || IsZero(right))
return MathExpressions.Zero;
if (IsOne(left))
return right;
if (IsOne(right))
return left;
break;

case ExpressionType.Divide:
if (IsZero(right))
throw new DivideByZeroException();
if (IsZero(left))
return MathExpressions.Zero;
if (IsOne(right))
return left;
break;

default:
throw new NotImplementedException();
}

return Expression.MakeBinary(node.NodeType, left, right);
}

protected override Expression VisitUnary(UnaryExpression node)
{
var operand = Visit(node.Operand);

ConstantExpression operandConstant = operand as ConstantExpression;
if (operandConstant != null && (operandConstant.Value is double))
{
double operandValue = (double)operandConstant.Value;

switch (node.NodeType)
{
case ExpressionType.Negate:
if (operandValue == 0.0)
return MathExpressions.Zero;

return Expression.Constant(-operandValue);

default:
throw new NotImplementedException();
}
}

switch (node.NodeType)
{
case ExpressionType.Negate:
if (operand.NodeType == ExpressionType.Negate)
{
return ((UnaryExpression)operand).Operand;
}

break;

default:
throw new NotImplementedException();
}

return Expression.MakeUnary(node.NodeType, operand, node.Type);
}

private static bool IsZero(Expression expression)
{
ConstantExpression constant = expression as ConstantExpression;
if (constant != null)
{
if (constant.Value.Equals(0.0))
return true;
}

return false;
}

private static bool IsOne(Expression expression)
{
ConstantExpression constant = expression as ConstantExpression;
if (constant != null)
{
if (constant.Value.Equals(1.0))
return true;
}

return false;
}
}
``````

## Formatting expressions for display with `ListPrintVisitor`

``````internal class ListPrintVisitor : ExpressionVisitor
{
protected override Expression VisitBinary(BinaryExpression node)
{
string op = null;

switch (node.NodeType)
{
op = "+";
break;
case ExpressionType.Subtract:
op = "-";
break;
case ExpressionType.Multiply:
op = "*";
break;
case ExpressionType.Divide:
op = "/";
break;
default:
throw new NotImplementedException();
}

var left = Visit(node.Left);
var right = Visit(node.Right);
string result = string.Format("({0} {1} {2})", op, ((ConstantExpression)left).Value, ((ConstantExpression)right).Value);
return Expression.Constant(result);
}

protected override Expression VisitConstant(ConstantExpression node)
{
if (node.Value is string)
return node;

return Expression.Constant(node.Value.ToString());
}

protected override Expression VisitParameter(ParameterExpression node)
{
return Expression.Constant(node.Name);
}
}
``````

## Testing the results

``````[TestMethod]
public void BasicSymbolicTest()
{
ParameterExpression x = Expression.Parameter(typeof(double), "x");
Expression linear = Expression.Add(Expression.Constant(3.0), x);
Assert.AreEqual("(+ 3 x)", Symbolic.ToString(linear));

Assert.AreEqual("(* (+ 3 x) (+ 2 x))", Symbolic.ToString(quadratic));

Expression expanded = Symbolic.Expand(quadratic);
Assert.AreEqual("(+ (+ (+ (* 3 2) (* 3 x)) (* x 2)) (* x x))", Symbolic.ToString(expanded));
Assert.AreEqual("(+ (+ (+ 6 (* 3 x)) (* x 2)) (* x x))", Symbolic.ToString(Symbolic.Simplify(expanded)));

Expression derivative = Symbolic.PartialDerivative(expanded, x);
Assert.AreEqual("(+ (+ (+ (+ (* 3 0) (* 0 2)) (+ (* 3 1) (* 0 x))) (+ (* x 0) (* 1 2))) (+ (* x 1) (* 1 x)))", Symbolic.ToString(derivative));

Expression simplified = Symbolic.Simplify(derivative);
Assert.AreEqual("(+ 5 (+ x x))", Symbolic.ToString(simplified));
}
``````
6
10/28/2009 3:28:22 AM

Prime Library

More Projects...