# symdiff.tcl --
#
#       Symbolic differentiation package for Tcl
#
# This package implements a single command, "math::symdiff::symdiff",
# which accepts a Tcl expression and a variable name, and if the expression
# is readily differentiable, returns a Tcl expression that evaluates the
# derivative.
#
# Copyright (c) 2005 by Kevin B. Kenny.  All rights reserved.
#
# See the file "license.terms" for information on usage and redistribution
# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
#
# RCS: @(#) $Id: symdiff.tcl,v not yet committed $


# This package requires the 'tclparser' from http://tclpro.sf.net/
# to analyze the expressions presented to it.

package require Tcl 8.4
package require parser 1.4
package provide math::symdiff 1.0

namespace eval math {}
namespace eval math::symdiff {
    namespace export symdiff
    namespace eval differentiate {}
}

# math::symdiff::symdiff --
#
#       Differentiate an expression with respect to a variable.
#
# Parameters:
#       expr -- expression to differentiate (Must be a Tcl expression,
#               without command substitution.)
#       var -- Name of the variable to differentiate the expression
#              with respect to.
#
# Results:
#       Returns a Tcl expression that evaluates the derivative.

proc math::symdiff::symdiff {expr var} {
    set parsetree [MakeParseTree $expr]
    return [ToInfix [differentiate::MakeDeriv $parsetree $var]]
}

# math::symdiff::MakeParseTree --
#
#       Build the parse tree for an expression
#
# Parameters:
#       expr -- Expression to parse
#
# Results:
#       Returns a parse tree.  The parse tree is a Tcl command
#       that can be evaluated in the math::symdiff::differentiate namespace
#       after the variable name is inserted at position 1.

proc math::symdiff::MakeParseTree {expr} {
    set parsetree [parse expr $expr [parse getrange $expr]]
    return [MakeParseTreeForNode $expr $parsetree]
}

# math::symdiff::MakeParseTreeForNode --
#
#       Converts a TclPro parse tree into one for the symbolic differentiator.
#
# Parameters:
#       expr -- Original expression
#       parseTree -- Parse tree from an expression, as returned by
#                    the 'parse' command in TclPro's 'tclparser'
#
# Results:
#       Returns a parse tree as for MakeParseTree.  Throws an error on anything
#       but a 'subexpression.'

proc math::symdiff::MakeParseTreeForNode {expr parsetree} {
    set kind [lindex $parsetree 0]
    switch -exact $kind {
        subexpr {
            return [MakeParseTreeForSubExpr $expr [lindex $parsetree 2]]
        }
        default {
            return -code error "symdiff doesn't handle $kind expressions"
        }
    }
}

# math::symdiff::MakeParseTreeForSubExpr --
#
#       Converts a TclPro parse tree into one for the symbolic differentiator.
#
# Parameters:
#       expr -- Original expression
#       parseTree -- Parse tree from one subexpression, as returned by
#                    the 'parse' command in TclPro's 'tclparser' Returned
#                    by [lindex $parseTree 2] if [lindex $parseTree 0] is
#                    'subexpr'
#
# Results:
#       Returns a parse tree as for MakeParseTree.  Throws an error on
#       command substitution (and anything else unrecognized).

proc math::symdiff::MakeParseTreeForSubExpr {expr components} {
    set kind [lindex $components 0 0]
    switch -exact $kind {
        operator {
            set tree [MakeOperator [parse getstring $expr \
                                       [lindex $components 0 1]]]
            foreach operand [lrange $components 1 end] {
                lappend tree [MakeParseTreeForNode $expr $operand]
            }
            return $tree
        }
        text {
            # a constant
            return [MakeConstant [parse getstring $expr \
                                      [lindex $components 0 1]]]
        }
        variable {
            return [MakeParseTreeForVariableExpr $expr \
                        [lindex $components 0 2]]
        }
        default {
            return -code error "symdiff doesn't handle $kind expressions"
        }
    }
}

# math::symdiff::MakeParseTreeForVariableExpr --
#
#       Converts a TclPro parse tree into one for the symbolic differentiator.
#
# Parameters:
#       expr -- Original expression
#       components -- Partial parse tree from TclPro's parser for the 
#                     variable reference.  Returned by [lrange $parsetree 2]
#                     if [lindex $parseTree 0] is 'variable'.
#
# Results:
#       Returns a parse tree as for MakeParseTree.  Throws an error on
#       array references.

proc math::symdiff::MakeParseTreeForVariableExpr {expr components} {
    if { [llength $components] == 1 } {
        set kind [lindex $components 0 0]
        if { $kind eq {text} } {
            set name [parse getstring $expr [lindex $components 0 1]]
            return [MakeVariable $name]
        } else {
            return -code error "symdiff expected a variable name, found $kind"
        }
    } else {
        return -code error "symdiff doesn't handle arrays"
    }
}

# math::symdiff::ToInfix --
#
#       Converts a parse tree to infix notation.
#
# Parameters:
#       tree - Parse tree to convert
#
# Results:
#       Returns the parse tree as a Tcl expression.

proc math::symdiff::ToInfix {tree} {
    set a [lindex $tree 0]
    set kind [lindex $a 0]
    switch -exact $kind {
        constant -
        text {
            set result [lindex $tree 1]
        }
        var {
            set result \$[lindex $tree 1]
        }
        operator {
            set name [lindex $a 1]
            if {([string is alnum $name] && $name ne {eq} && $name ne {ne})
                || [llength $tree] == 2} {
                set result $name
                append result \(
                set sep ""
                foreach arg [lrange $tree 1 end] {
                    append result $sep [ToInfix $arg]
                    set sep ", "
                }
                append result \)
            } elseif {[llength $tree] == 3} {
                set result \(
                append result [ToInfix [lindex $tree 1]]
                append result " " $name " "
                append result [ToInfix [lindex $tree 2]]
                append result \)
            } else {
                error "symdiff encountered a malformed parse, can't happen"
            }
        }
        default {
            error "symdiff can't synthesize a $kind expression"
        }
    }
    return $result
}

# math::symdiff::differentiate::MakeDeriv --
#
#       Differentiates a Tcl expression represented as a parse tree.
#
# Parameters:
#       tree -- Parse tree from MakeParseTreeForExpr
#       var -- Variable to differentiate with respect to
#
# Results:
#       Returns the parse tree of the derivative.

proc math::symdiff::differentiate::MakeDeriv {tree var} {
    return [eval [linsert $tree 1 $var]]
}

# math::symdiff::differentiate::ChainRule --
#
#       Applies the Chain Rule to evaluate the derivative of a unary 
#       function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       derivMaker -- Command prefix for differentiating the function.
#       u -- Function argument.
#
# Results:
#       Returns a parse tree representing the derivative of f($u).
#
# ChainRule differentiates $u with respect to $var by calling MakeDeriv,
# makes the derivative of f($u) with respect to $u by calling derivMaker
# passing $u as a parameter, and then returns a parse tree representing
# the product of the results.

proc math::symdiff::differentiate::ChainRule {var derivMaker u} {
    lappend derivMaker $u
    set result [MakeProd [MakeDeriv $u $var] [eval $derivMaker]]
}

# math::symdiff::differentiate::constant --
#
#       Differentiate a constant.
#
# Parameters:
#       var -- Variable to differentiate with respect to - unused
#       constant -- Constant expression to differentiate - ignored
#
# Results:
#       Returns a parse tree of the derivative, which is, of course, the
#       constant zero.

proc math::symdiff::differentiate::constant {var constant} {
    return [MakeConstant 0.0]
}

# math::symdiff::differentiate::var --
#
#       Differentiate a variable expression.
#
# Parameters:
#       var - Variable with which to differentiate.
#       exprVar - Expression being differentiated, which is a single
#                 variable.
#
# Results:
#       Returns a parse tree of the derivative.
#
# The derivative is the constant unity if the variables are the same
# and the constant zero if they are different.

proc math::symdiff::differentiate::var {var exprVar} {
    if { $exprVar eq $var } {
        return [MakeConstant 1.0]
    } else {
        return [MakeConstant 0.0]
    }
}

# math::symdiff::differentiate::operator + --
#
#       Forms the derivative of a sum.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       args -- One or two arguments giving augend and addend. If only
#               one argument is supplied, this is unary +.
#
# Results:
#       Returns a parse tree representing the derivative.
#
# Of course, the derivative of a sum is the sum of the derivatives.

proc {math::symdiff::differentiate::operator +} {var args} {
    if { [llength $args] == 1 } {
        set u [lindex $args 0]
        set result [eval [linsert $u 1 $var]]
    } elseif { [llength $args] == 2 } {
        foreach {u v} $args break
        set du [eval [linsert $u 1 $var]]
        set dv [eval [linsert $v 1 $var]]
        set result [MakeSum $du $dv]
    } else {
        error "symdiff encountered a malformed parse, can't happen"
    }
    return $result
}

# math::symdiff::differentiate::operator - --
#
#       Forms the derivative of a difference.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       args -- One or two arguments giving minuend and subtrahend. If only
#               one argument is supplied, this is unary -.
#
# Results:
#       Returns a parse tree representing the derivative.
#
# Of course, the derivative of a sum is the sum of the derivatives.

proc {math::symdiff::differentiate::operator -} {var args} {
    if { [llength $args] == 1 } {
        set u [lindex $args 0]
        set du [eval [linsert $u 1 $var]]
        set result [MakeUnaryMinus $du]
    } elseif { [llength $args] == 2 } {
        foreach {u v} $args break
        set du [eval [linsert $u 1 $var]]
        set dv [eval [linsert $v 1 $var]]
        set result [MakeDifference $du $dv]
    } else {
        error "symdiff encounered a malformed parse, can't happen"
    }
    return $result
}

# math::symdiff::differentiate::operator * --
#
#       Forms the derivative of a product.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u, v -- Multiplicand and multiplier. 
#
# Results:
#       Returns a parse tree representing the derivative.
#
# The familiar freshman calculus product rule.

proc {math::symdiff::differentiate::operator *} {var u v} {
    set du [eval [linsert $u 1 $var]]
    set dv [eval [linsert $v 1 $var]]
    set result [MakeSum [MakeProd $dv $u] [MakeProd $du $v]]
    return $result
}

# math::symdiff::differentiate::operator / --
#
#       Forms the derivative of a quotient.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u, v -- Dividend and divisor. 
#
# Results:
#       Returns a parse tree representing the derivative.
#
# The familiar freshman calculus quotient rule.

proc {math::symdiff::differentiate::operator /} {var u v} {
    set du [eval [linsert $u 1 $var]]
    set dv [eval [linsert $v 1 $var]]
    set result [MakeQuotient \
                    [MakeDifference \
                         $du \
                         [MakeQuotient \
                              [MakeProd $dv $u] \
                              $v]] \
                    $v]
    return $result
}

# math::symdiff::differentiate::operator acos --
#
#       Differentiates the 'acos' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the acos() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain Rule: D(acos(u))=-D(u)/sqrt(1 - u*u)
# (Might it be better to factor 1-u*u into (1+u)(1-u)? Less likely to be
# catastrophic cancellation if u is near 1?)

proc {math::symdiff::differentiate::operator acos} {var u} {
    set du [eval [linsert $u 1 $var]]
    set result [MakeQuotient [MakeUnaryMinus $du] \
                    [MakeFunCall sqrt \
                         [MakeDifference [MakeConstant 1.0] \
                              [MakeProd $u $u]]]]
    return $result
}

# math::symdiff::differentiate::operator asin --
#
#       Differentiates the 'asin' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the asin() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain Rule: D(asin(u))=D(u)/sqrt(1 - u*u)
# (Might it be better to factor 1-u*u into (1+u)(1-u)? Less likely to be
# catastrophic cancellation if u is near 1?)

proc {math::symdiff::differentiate::operator asin} {var u} {
    set du [eval [linsert $u 1 $var]]
    set result [MakeQuotient $du \
                    [MakeFunCall sqrt \
                         [MakeDifference [MakeConstant 1.0] \
                              [MakeProd $u $u]]]]
    return $result
}

# math::symdiff::differentiate::operator atan --
#
#       Differentiates the 'atan' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the atan() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain Rule: D(atan(u))=D(u)/(1 + $u*$u)

proc {math::symdiff::differentiate::operator atan} {var u} {
    set du [eval [linsert $u 1 $var]]
    set result [MakeQuotient $du \
                    [MakeSum [MakeConstant 1.0] \
                         [MakeProd $u $u]]]
}

# math::symdiff::differentiate::operator atan2 --
#
#       Differentiates the 'atan2' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       f, g -- Arguments to the atan() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain and Quotient Rules: 
#       D(atan2(f, g)) = (D(f)*g - D(g)*f)/(f*f + g*g)

proc {math::symdiff::differentiate::operator atan2} {var f g} {
    set df [eval [linsert $f 1 $var]]
    set dg [eval [linsert $g 1 $var]]
    return [MakeQuotient \
                [MakeDifference \
                     [MakeProd $df $g] \
                     [MakeProd $f $dg]] \
                [MakeSum \
                     [MakeProd $f $f] \
                     [MakeProd $g $g]]]
}

# math::symdiff::differentiate::operator cos --
#
#       Differentiates the 'cos' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the cos() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain Rule: D(cos(u))=-sin(u)*D(u)

proc {math::symdiff::differentiate::operator cos} {var u} {
    return [ChainRule $var MakeMinusSin $u]
}
proc math::symdiff::differentiate::MakeMinusSin {operand} {
    return [MakeUnaryMinus [MakeFunCall sin $operand]]
}

# math::symdiff::differentiate::operator cosh --
#
#       Differentiates the 'cosh' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the cosh() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain Rule: D(cosh(u))=sinh(u)*D(u)

proc {math::symdiff::differentiate::operator cosh} {var u} {
    set result [ChainRule $var [list MakeFunCall sinh] $u]
    return $result
}

# math::symdiff::differentiate::operator exp --
#
#       Differentiate the exponential function
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument of the exponential function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Uses the Chain Rule D(exp(u)) = exp(u)*D(u).

proc {math::symdiff::differentiate::operator exp} {var u} {
    set result [ChainRule $var [list MakeFunCall exp] $u]
    return $result
}

# math::symdiff::differentiate::operator hypot --
#
#       Differentiate the 'hypot' function
#
# Parameters:
#       var - Variable to differentiate with respect to.
#       f, g - Arguments to the 'hypot' function
#
# Results:
#       Returns a parse tree of the derivative
#
# Uses a number of algebraic simplifications to arrive at:
#       D(hypot(f,g)) = (f*D(f)+g*D(g))/hypot(f,g)

proc {math::symdiff::differentiate::operator hypot} {var f g} {
    set df [eval [linsert $f 1 $var]]
    set dg [eval [linsert $g 1 $var]]
    return [MakeQuotient \
                [MakeSum \
                     [MakeProd $df $f] \
                     [MakeProd $dg $g]] \
                [MakeFunCall hypot $f $g]]
}

# math::symdiff::differentiate::operator log --
#
#       Differentiates a logarithm.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the log() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# D(log(u))==D(u)/u

proc {math::symdiff::differentiate::operator log} {var u} {
    set du [eval [linsert $u 1 $var]]
    set result [MakeQuotient $du $u]
    return $result
}

# math::symdiff::differentiate::operator log10 --
#
#       Differentiates a common logarithm.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the log10() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# D(log(u))==D(u)/(u * log(10))

proc {math::symdiff::differentiate::operator log10} {var u} {
    set du [eval [linsert $u 1 $var]]
    set result [MakeQuotient $du \
                    [MakeProd [MakeConstant [expr log(10.)]] $u]]
    return $result
}

# math::symdiff::differentiate::operator pow --
#
#       Differentiate an exponential.
#
# Parameters:
#       var -- Variable to differentiate with respect to
#       f, g -- Base and exponent
#
# Results:
#       Returns the parse tree of the derivative.
#
# Handles the special case where g is constant as
#    D(f**g) == g*f**(g-1)*D(f)
# Otherwise, uses the general power formula
#    D(f**g) == (f**g) * (((D(f)*g)/f) + (D(g)*log(f)))

proc {math::symdiff::differentiate::operator pow} {var f g} {
    set df [eval [linsert $f 1 $var]]
    if { [IsConstant $g] } {
        set gm1 [MakeConstant [expr {[ConstantValue $g] - 1}]]
        set result [MakeProd $df [MakeProd $g [MakePower $f $gm1]]]
        
    } else {
        set dg [eval [linsert $g 1 $var]]
        set result [MakeProd [MakePower $f $g] \
                        [MakeSum \
                             [MakeQuotient [MakeProd $df $g] $f] \
                             [MakeProd $dg [MakeFunCall log $f]]]]
    }
    return $result
}

# math::symdiff::differentiate::operator sin --
#
#       Differentiates the 'sin' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the sin() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain Rule: D(sin(u))=cos(u)*D(u)

proc {math::symdiff::differentiate::operator sin} {var u} {
    set result [ChainRule $var [list MakeFunCall cos] $u]
    return $result
}

# math::symdiff::differentiate::operator sinh --
#
#       Differentiates the 'sinh' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the sinh() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain Rule: D(sin(u))=cosh(u)*D(u)

proc {math::symdiff::differentiate::operator sinh} {var u} {
    set result [ChainRule $var [list MakeFunCall cosh] $u]
    return $result
}

# math::symdiff::differentiate::operator sqrt --
#
#       Differentiate the 'sqrt' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to
#       u -- Parameter of 'sqrt' as a parse tree.
#
# Results:
#       Returns a parse tree representing the derivative.
#
# D(sqrt(u))==D(u)/(2*sqrt(u))

proc {math::symdiff::differentiate::operator sqrt} {var u} {
    set du [eval [linsert $u 1 $var]]
    set result [MakeQuotient $du [MakeProd [MakeConstant 2.0] \
                                      [MakeFunCall sqrt $u]]]
    return $result
}

# math::symdiff::differentiate::operator tan --
#
#       Differentiates the 'tan' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the tan() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain Rule: D(tan(u))=D(u)/(cos(u)*cos(u))

proc {math::symdiff::differentiate::operator tan} {var u} {
    set du [eval [linsert $u 1 $var]]
    set cosu [MakeFunCall cos $u]
    return [MakeQuotient $du [MakeProd $cosu $cosu]]
}

# math::symdiff::differentiate::operator tanh --
#
#       Differentiates the 'tanh' function.
#
# Parameters:
#       var -- Variable to differentiate with respect to.
#       u -- Argument to the tanh() function.
#
# Results:
#       Returns a parse tree of the derivative.
#
# Applies the Chain Rule: D(tanh(u))=D(u)/(cosh(u)*cosh(u))

proc {math::symdiff::differentiate::operator tanh} {var u} {
    set du [eval [linsert $u 1 $var]]
    set coshu [MakeFunCall cosh $u]
    return [MakeQuotient $du [MakeProd $coshu $coshu]]
}

# math::symdiff::MakeFunCall --
#
#       Makes a parse tree for a function call
#
# Parameters:
#       fun -- Name of the function to call
#       args -- Arguments to the function, expressed as parse trees
#
# Results:
#       Returns a parse tree for the result of calling the function.
#
# Performs the peephole optimization of replacing a function with
# constant parameters with its value.

proc math::symdiff::MakeFunCall {fun args} {
    set constant 1
    set exp $fun
    append exp \(
    set sep ""
    foreach a $args {
        if {[IsConstant $a]} {
            append exp $sep [ConstantValue $a]
            set sep ","
        } else {
            set constant 0
            break
        }
    }
    if { $constant } {
        append exp \)
        return [MakeConstant [expr $exp]]
    }
    set result [MakeOperator $fun]
    foreach arg $args {
        lappend result $arg
    }
    return $result
}

# math::symdiff::MakeSum --
#
#       Makes the parse tree for a sum.
#
# Parameters:
#       left, right -- Parse trees for augend and addend
#
# Results:
#       Returns the parse tree for the sum.
#
# Performs the following peephole optimizations:
# (1) a + (-b) = a - b
# (2) (-a) + b = b - a
# (3) 0 + a = a
# (4) a + 0 = a
# (5) The sum of two constants may be reduced to a constant

proc math::symdiff::MakeSum {left right} {
    if { [IsUnaryMinus $right] } {
        return [MakeDifference $left [UnaryMinusArg $right]]
    }
    if { [IsUnaryMinus $left] } {
        return [MakeDifference $right [UnaryMinusArg $left]]
    }
    if { [IsConstant $left] } {
        set v [ConstantValue $left]
        if { $v == 0 } {
            return $right
        } elseif { [IsConstant $right] } {
            return [MakeConstant [expr {[ConstantValue $left]
                                        + [ConstantValue $right]}]]
        }
    } elseif { [IsConstant $right] } {
        set v [ConstantValue $right]
        if { $v == 0 } {
            return $left
        }
    }
    set result [MakeOperator +]
    lappend result $left $right
    return $result
}

# math::symdiff::MakeDifference --
#
#       Makes the parse tree for a difference
#
# Parameters:
#       left, right -- Minuend and subtrahend, expressed as parse trees
#
# Results:
#       Returns a parse tree expressing the difference
#
# Performs the following peephole optimizations:
# (1) a - (-b) = a + b
# (2) -a - b = -(a + b)
# (3) 0 - b = -b
# (4) a - 0 = a
# (5) The difference of any two constants can be reduced to a constant.

proc math::symdiff::MakeDifference {left right} {
    if { [IsUnaryMinus $right] } {
        return [MakeSum $left [UnaryMinusArg $right]]
    }
    if { [IsUnaryMinus $left] } {
        return [MakeUnaryMinus [MakeSum [UnaryMinusArg $left] $right]]
    }
    if { [IsConstant $left] } {
        set v [ConstantValue $left]
        if { $v == 0 } {
            return [MakeUnaryMinus $right]
        } elseif { [IsConstant $right] } {
            return [MakeConstant [expr {[ConstantValue $left]
                                        - [ConstantValue $right]}]]
        }
    } elseif { [IsConstant $right] } {
        set v [ConstantValue $right]
        if { $v == 0 } {
            return $left
        }
    }
    set result [MakeOperator -]
    lappend result $left $right
    return $result
}

# math::symdiff::MakeProd --
#
#       Constructs the parse tree for a product, left*right.
#
# Parameters:
#       left, right - Multiplicand and multiplier
#
# Results:
#       Returns the parse tree for the result.
#
# Performs the following peephole optimizations.
# (1) If either operand is a unary minus, it is hoisted out of the
#     expression.
# (2) If either operand is the constant 0, the result is the constant 0
# (3) If either operand is the constant 1, the result is the other operand.
# (4) If either operand is the constant -1, the result is unary minus
#     applied to the other operand
# (5) If both operands are constant, the result is a constant containing
#     their product.

proc math::symdiff::MakeProd {left right} {
    if { [IsUnaryMinus $left] } {
        return [MakeUnaryMinus [MakeProd [UnaryMinusArg $left] $right]]
    }
    if { [IsUnaryMinus $right] } {
        return [MakeUnaryMinus [MakeProd $left [UnaryMinusArg $right]]]
    }
    if { [IsConstant $left] } {
        set v [ConstantValue $left]
        if { $v == 0 } {
            return [MakeConstant 0.0]
        } elseif { $v == 1 } {
            return $right
        } elseif { $v == -1 } {
            return [MakeUnaryMinus $right]
        } elseif { [IsConstant $right] } {
            return [MakeConstant [expr {[ConstantValue $left]
                                        * [ConstantValue $right]}]]
        }
    } elseif { [IsConstant $right] } {
        set v [ConstantValue $right]
        if { $v == 0 } {
            return [MakeConstant 0.0]
        } elseif { $v == 1 } {
            return $left
        } elseif { $v == -1 } {
            return [MakeUnaryMinus $left]
        }
    }
    set result [MakeOperator *]
    lappend result $left $right
    return $result
}

# math::symdiff::MakeQuotient --
#
#       Makes a parse tree for a quotient, n/d
#
# Parameters:
#       n, d - Parse trees for numerator and denominator
#
# Results:
#       Returns the parse tree for the quotient.
#
# Performs peephole optimizations:
# (1) If either operand is a unary minus, it is hoisted out.
# (2) If the numerator is the constant 0, the result is the constant 0.
# (3) If the demominator is the constant 1, the result is the numerator
# (4) If the denominator is the constant -1, the result is the unary
#     negation of the numerator.
# (5) If both numerator and denominator are constant, the result is
#     a constant representing their quotient.

proc math::symdiff::MakeQuotient {n d} {
    if { [IsUnaryMinus $n] } {
        return [MakeUnaryMinus [MakeQuotient [UnaryMinusArg $n] $d]]
    }
    if { [IsUnaryMinus $d] } {
        return [MakeUnaryMinus [MakeQuotient $n [UnaryMinusArg $d]]]
    }
    if { [IsConstant $n] } {
        set v [ConstantValue $n]
        if { $v == 0 } {
            return [MakeConstant 0.0]
        } elseif { [IsConstant $d] } {
            return [MakeConstant [expr {[ConstantValue $n]
                                        * [ConstantValue $d]}]]
        }
    } elseif { [IsConstant $d] } {
        set v [ConstantValue $d]
        if { $v == 0 } {
            return -code error "requested expression will result in division by zero at run time"
        } elseif { $v == 1 } {
            return $n
        } elseif { $v == -1 } {
            return [MakeUnaryMinus $n]
        }
    }
    set result [MakeOperator /]
    lappend result $n $d
    return $result
}

# math::symdiff::MakePower --
#
#       Make a parse tree for an exponentiation operation
#
# Parameters:
#       a -- Base, expressed as a parse tree
#       b -- Exponent, expressed as a parse tree
#
# Results:
#       Returns a parse tree for the expression
#
# Performs peephole optimizations:
# (1) The constant zero raised to any non-zero power is 0
# (2) The constant 1 raised to any power is 1
# (3) Any non-zero quantity raised to the zero power is 1
# (4) Any non-zero quantity raised to the first power is the base itself.
# (5) MakeFunCall will optimize any other case of a constant raised
#     to a constant power.

proc math::symdiff::MakePower {a b} {
    if { [IsConstant $a] } {
        if { [ConstantValue $a] == 0 } {
            if { [IsConstant $b] && [ConstantValue $b] == 0 } {
                error "requested expression will result in zero to zero power at run time"
            }
            return [MakeConstant 0.0]
        } elseif { [ConstantValue $a] == 1 } {
            return [MakeConstant 1.0]
        }
    }
    if { [IsConstant $b] } {
        if { [ConstantValue $b] == 0 } {
            return [MakeConstant 1.0]
        } elseif { [ConstantValue $b] == 1 } {
            return $a
        }
    }
    return [MakeFunCall pow $a $b]
}

# math::symdiff::MakeUnaryMinus --
#
#       Makes the parse tree for a unary negation.
#
# Parameters:
#       operand -- Parse tree for the operand
#
# Results:
#       Returns the parse tree for the expression
#
# Performs the following peephole optimizations:
# (1) -(-$a) = $a
# (2) The unary negation of a constant is another constant

proc math::symdiff::MakeUnaryMinus {operand} {
    if { [IsUnaryMinus $operand] } {
        return [UnaryMinusArg $operand]
    }
    if { [IsConstant $operand] } {
        return [MakeConstant [expr {-[ConstantValue $operand]}]]
    } else {
        return [list [list operator -] $operand]
    }
}

# math::symdiff::IsUnaryMinus --
#
#       Determines whether a parse tree represents a unary negation
#
# Parameters:
#       x - Parse tree to examine
#
# Results:
#       Returns 1 if the parse tree represents a unary minus, 0 otherwise

proc math::symdiff::IsUnaryMinus {x} {
    return [expr {[llength $x] == 2
                  && [lindex $x 0] eq [list operator -]}]
}

# math::symdiff::UnaryMinusArg --
#
#       Extracts the argument from a unary negation.
#
# Parameters:
#       x - Parse tree to examine, known to represent a unary negation
#
# Results:
#       Returns a parse tree representing the operand.

proc math::symdiff::UnaryMinusArg {x} {
    return [lindex $x 1]
}

# math::symdiff::MakeOperator --
#
#       Makes a partial parse tree for an operator
#
# Parameters:
#       op -- Name of the operator
#
# Results:
#       Returns the resulting parse tree.
#
# The caller may use [lappend] to place any needed operands

proc math::symdiff::MakeOperator {op} {
    if { $op eq {?} } {
        return -code error "symdiff can't differentiate the ternary ?: operator"
    } elseif { [llength [info commands [list differentiate::operator $op]]] != 0 } {
        return [list [list operator $op]]
    } elseif { [string is alnum $op] && ($op ne {eq}) && ($op ne {ne}) } {
        return -code error "symdiff can't differentiate the \"$op\" function"
    } else {
        return -code error "symdiff can't differentiate the \"$op\" operator"
    }
}

# math::symdiff::MakeVariable --
#
#       Makes a partial parse tree for a single variable
#
# Parameters:
#       name -- Name of the variable
#
# Results:
#       Returns a partial parse tree giving the variable

proc math::symdiff::MakeVariable {name} {
    return [list var $name]
}

# math::symdiff::MakeConstant --
#
#       Make the parse tree for a constant.
#
# Parameters:
#       value -- The constant's value
#
# Results:
#       Returns a parse tree.

proc math::symdiff::MakeConstant {value} {
    return [list constant $value]
}

# math::symdiff::IsConstant --
#
#       Test if an expression represented by a parse tree is a constant.
#
# Parameters:
#       Item - Parse tree to test
#
# Results:
#       Returns 1 for a constant, 0 for anything else

proc math::symdiff::IsConstant {item} {
    return [expr {[lindex $item 0] eq {constant}}]
}

# math::symdiff::ConstantValue --
#
#       Recovers a constant value from the parse tree representing a constant
#       expression.
#
# Parameters:
#       item -- Parse tree known to be a constant.
#
# Results:
#       Returns the constant value.

proc math::symdiff::ConstantValue {item} {
    return [lindex $item 1]
}

# Define the parse tree fabrication routines in the 'differentiate'
# namespace as well as the 'symdiff' namespace, without exporting them
# from the package.

interp alias {} math::symdiff::differentiate::IsConstant \
    {} math::symdiff::IsConstant
interp alias {} math::symdiff::differentiate::ConstantValue \
    {} math::symdiff::ConstantValue
interp alias {} math::symdiff::differentiate::MakeConstant \
    {} math::symdiff::MakeConstant
interp alias {} math::symdiff::differentiate::MakeDifference \
    {} math::symdiff::MakeDifference
interp alias {} math::symdiff::differentiate::MakeFunCall \
    {} math::symdiff::MakeFunCall
interp alias {} math::symdiff::differentiate::MakePower \
    {} math::symdiff::MakePower
interp alias {} math::symdiff::differentiate::MakeProd \
    {} math::symdiff::MakeProd
interp alias {} math::symdiff::differentiate::MakeQuotient \
    {} math::symdiff::MakeQuotient
interp alias {} math::symdiff::differentiate::MakeSum \
    {} math::symdiff::MakeSum
interp alias {} math::symdiff::differentiate::MakeUnaryMinus \
    {} math::symdiff::MakeUnaryMinus
interp alias {} math::symdiff::differentiate::ExtractExpression \
    {} math::symdiff::ExtractExpression
