多项式的评价生成的方法

问题描述

我试图想出一个优雅的方式来处理一些生成多项式。以下是我们将重点关注（独家）为这一问题的情况：

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. 的为了的是在生成参数一个的 N 的阶多项式，其中n：=秩序+ 1


2. 的 I 的取值范围为0的整数参数.. ñ


3. 多项式具有x_j，其中j = 1..N和j≠我（应该是在这一点上需要计算器一项新功能，清晰的零或它的存在，我不知道吧）


4. 多项式的计算结果为1的x_i。

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.

由于这种特殊的代码示例产生X_1 ... x_n，我将解释他们是如何在代码中找到。这些点间隔均匀的 x_j = j的* elementSize /订单分开，其中 N =订单+ 1 。

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.

我生成一个 Func键<双层，双方式> 来评估这个polynomial¹

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;
        factors.Add(Expression.Subtract(Expression.Constant(p), xp));
    }

    // 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();
}

问题：我还需要评估ψ '=dψ/ DX。要做到这一点，我可以重写ψ=×规模（X_0 - X）（X_1 - X）×..×（x_n - X）/（x_i - x）的形式ψ=α_n×X ^ N +α_ntimes; X ^（N-1）+ ... +α_1×X +α_0。这使ψ'= N×α_n×X ^（N-1）+（N-1）×α_n×X ^（N-2）+ ... + 1×α_1

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.

有关计算的原因，我们可以重写写ψ= X×（X×（X×（不调用 Math.Pow 最后的答案.. ） - β_2） - β_1） - β_0

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. 展开一个因素表达式包含常量表达式和 ParameterExpression 树叶和基本的数学运算（落得无论是 BinaryExpression 与的NodeType 设置为操作） - 结果在这里可以包括 InvocationExpression 元素的的MethodInfo Math.Pow ，我们将在整个一个特殊的方式来处理。


2. 然后我把导相对于一些指定的 ParameterExpression 。在条款的结果，其中右侧参数 Math.Pow 的调用是常数2是由常量表达式（2）乘什么左侧（。当调用的右手边就变成了常量表达式与值 1 ，我们与更换调用只是 ParameterExpression 本身。



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.

¹在未来，我想的方法采取 ParameterExpression 并返回的表达式，其基于该参数。这样我可以聚合生成的函数。我还没有。
²在未来，我希望能发布一个通用库与LINQ表达式作为象征性的数学工作。

¹ 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.

推荐答案

我使用的键入.NET 4它并不完美，但它看起来像一个可行的解决方案的基础。

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.

• 符号是一个公共静态类暴露的方法，如扩展，简化和 PartialDerivative


• ExpandVisitor 是一个内部的辅助类型扩展表达式


• SimplifyVisitor 是一个简化的表达式内部helper类


• DerivativeVisitor 是需要表达的微分内部helper类


• ListPrintVisitor 是内部，一个表达式转换为一个前缀表示了Lisp语言的语法

型辅助

• 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

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();
    }
}





扩大表达与 ExpandVisitor

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)
            {
            case ExpressionType.Add:
                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;
    }
}

隔空衍生 DerivativeVisitor

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.Add:
        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
            {
                return Expression.Multiply(node, Expression.Add(
                    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;
    }
}

简化表达与 SimplifyVisitor

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)
            {
            case ExpressionType.Add:
                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)
        {
        case ExpressionType.Add:
            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;
    }
}

格式表达式显示与 ListPrintVisitor

Formatting expressions for display with ListPrintVisitor

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

        switch (node.NodeType)
        {
        case ExpressionType.Add:
            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));

    Expression quadratic = Expression.Multiply(linear, Expression.Add(Expression.Constant(2.0), x));
    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));
}

这篇关于多项式的评价生成的方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！