1# symdiff.tcl -- 2# 3# Symbolic differentiation package for Tcl 4# 5# This package implements a command, "math::calculus::symdiff::symdiff", 6# which accepts a Tcl expression and a variable name, and if the expression 7# is readily differentiable, returns a Tcl expression that evaluates the 8# derivative. 9# 10# Copyright (c) 2005, 2010 by Kevin B. Kenny. All rights reserved. 11# 12# See the file "license.terms" for information on usage and redistribution 13# of this file, and for a DISCLAIMER OF ALL WARRANTIES. 14# 15# RCS: @(#) $Id: symdiff.tcl,v 1.2 2011/01/13 02:49:53 andreas_kupries Exp $ 16 17 18# This package requires the 'tclparser' from http://tclpro.sf.net/ 19# to analyze the expressions presented to it. 20 21package require Tcl 8.4 22package require grammar::aycock 1.0 23package provide math::calculus::symdiff 1.0.1 24 25namespace eval math {} 26namespace eval math::calculus {} 27namespace eval math::calculus::symdiff { 28 namespace export jacobian symdiff 29 namespace eval differentiate {} 30} 31 32# math::calculus::symdiff::jacobian -- 33# 34# Differentiate a set of expressions with respect to a set of 35# model variables 36# 37# Parameters: 38# model -- A list of alternating {variable name} {expr} 39# 40# Results: 41# Returns a list of lists. The ith sublist is the gradient vector 42# of the ith expr in the model; that is, the jth element of 43# the ith sublist is the derivative of the ith expr with respect 44# to the jth variable. 45# 46# Returns an error if any expression cannot be differentiated with 47# respect to any of the elements of the list, or if the list has 48# no elements or an odd number of elements. 49 50proc math::calculus::symdiff::jacobian {list} { 51 set l [llength $list] 52 if {$l == 0 || $l%2 != 0} { 53 return -code error "list of variables and expressions must have an odd number of elements" 54 } 55 set J {} 56 foreach {- expr} $list { 57 set gradient {} 58 foreach {var -} $list { 59 lappend gradient [symdiff $expr $var] 60 } 61 lappend J $gradient 62 } 63 return $J 64} 65 66# math::calculus::symdiff::symdiff -- 67# 68# Differentiate an expression with respect to a variable. 69# 70# Parameters: 71# expr -- expression to differentiate (Must be a Tcl expression, 72# without command substitution.) 73# var -- Name of the variable to differentiate the expression 74# with respect to. 75# 76# Results: 77# Returns a Tcl expression that evaluates the derivative. 78 79proc math::calculus::symdiff::symdiff {expr var} { 80 variable parser 81 set parsetree [$parser parse {*}[Lexer $expr] [namespace current]] 82 return [ToInfix [differentiate::MakeDeriv $parsetree $var]] 83} 84 85# math::calculus::symdiff::Parser -- 86# 87# Parser for the mathematical expressions that this package can 88# differentiate. 89 90namespace eval math::calculus::symdiff { 91 variable parser [grammar::aycock::parser { 92 expression ::= expression addop term { 93 set result [${clientData}::MakeOperator [lindex $_ 1]] 94 lappend result [lindex $_ 0] [lindex $_ 2] 95 } 96 expression ::= term { 97 lindex $_ 0 98 } 99 100 addop ::= + { 101 lindex $_ 0 102 } 103 addop ::= - { 104 lindex $_ 0 105 } 106 107 term ::= term mulop factor { 108 set result [${clientData}::MakeOperator [lindex $_ 1]] 109 lappend result [lindex $_ 0] [lindex $_ 2] 110 } 111 term ::= factor { 112 lindex $_ 0 113 } 114 mulop ::= * { 115 lindex $_ 0 116 } 117 mulop ::= / { 118 lindex $_ 0 119 } 120 121 factor ::= addop factor { 122 set result [${clientData}::MakeOperator [lindex $_ 0]] 123 lappend result [lindex $_ 1] 124 } 125 factor ::= expon { 126 lindex $_ 0 127 } 128 129 expon ::= primary ** expon { 130 set result [${clientData}::MakeOperator [lindex $_ 1]] 131 lappend result [lindex $_ 0] [lindex $_ 2] 132 } 133 expon ::= primary { 134 lindex $_ 0 135 } 136 137 primary ::= {$} bareword { 138 ${clientData}::MakeVariable [lindex $_ 1] 139 } 140 primary ::= number { 141 ${clientData}::MakeConstant [lindex $_ 0] 142 } 143 primary ::= bareword ( arglist ) { 144 set result [${clientData}::MakeOperator [lindex $_ 0]] 145 lappend result {*}[lindex $_ 2] 146 } 147 primary ::= ( expression ) { 148 lindex $_ 1 149 } 150 151 arglist ::= expression { 152 set _ 153 } 154 arglist ::= arglist , expression { 155 linsert [lindex $_ 0] end [lindex $_ 2] 156 } 157 }] 158} 159 160# math::calculus::symdiff::Lexer -- 161# 162# Lexer for the arithmetic expressions that the 'symdiff' package 163# can differentiate. 164# 165# Results: 166# Returns a two element list. The first element is a list of the 167# lexical values of the tokens that were found in the expression; 168# the second is a list of the semantic values of the tokens. The 169# two sublists are the same length. 170 171proc math::calculus::symdiff::Lexer {expression} { 172 set start 0 173 set tokens {} 174 set values {} 175 while {$expression ne {}} { 176 if {[regexp {^\*\*(.*)} $expression -> rest]} { 177 178 # Exponentiation 179 180 lappend tokens ** 181 lappend values ** 182 } elseif {[regexp {^([-+/*$(),])(.*)} $expression -> token rest]} { 183 184 # Single-character operators 185 186 lappend tokens $token 187 lappend values $token 188 } elseif {[regexp {^([[:alpha:]][[:alnum:]_]*)(.*)} \ 189 $expression -> token rest]} { 190 191 # Variable and function names 192 193 lappend tokens bareword 194 lappend values $token 195 } elseif {[regexp -nocase -expanded { 196 ^((?: 197 (?: [[:digit:]]+ (?:[.][[:digit:]]*)? ) 198 | (?: [.][[:digit:]]+ ) ) 199 (?: e [-+]? [[:digit:]]+ )? ) 200 (.*) 201 }\ 202 $expression -> token rest]} { 203 204 # Numbers 205 206 lappend tokens number 207 lappend values $token 208 } elseif {[regexp {^[[:space:]]+(.*)} $expression -> rest]} { 209 210 # Whitespace 211 212 } else { 213 214 # Anything else is an error 215 216 return -code error \ 217 -errorcode [list MATH SYMDIFF INVCHAR \ 218 [string index $expression 0]] \ 219 [list invalid character [string index $expression 0]] \ 220 } 221 set expression $rest 222 } 223 return [list $tokens $values] 224} 225 226# math::calculus::symdiff::ToInfix -- 227# 228# Converts a parse tree to infix notation. 229# 230# Parameters: 231# tree - Parse tree to convert 232# 233# Results: 234# Returns the parse tree as a Tcl expression. 235 236proc math::calculus::symdiff::ToInfix {tree} { 237 set a [lindex $tree 0] 238 set kind [lindex $a 0] 239 switch -exact $kind { 240 constant - 241 text { 242 set result [lindex $tree 1] 243 } 244 var { 245 set result \$[lindex $tree 1] 246 } 247 operator { 248 set name [lindex $a 1] 249 if {([string is alnum $name] && $name ne {eq} && $name ne {ne}) 250 || [llength $tree] == 2} { 251 set result $name 252 append result \( 253 set sep "" 254 foreach arg [lrange $tree 1 end] { 255 append result $sep [ToInfix $arg] 256 set sep ", " 257 } 258 append result \) 259 } elseif {[llength $tree] == 3} { 260 set result \( 261 append result [ToInfix [lindex $tree 1]] 262 append result " " $name " " 263 append result [ToInfix [lindex $tree 2]] 264 append result \) 265 } else { 266 error "symdiff encountered a malformed parse, can't happen" 267 } 268 } 269 default { 270 error "symdiff can't synthesize a $kind expression" 271 } 272 } 273 return $result 274} 275 276# math::calculus::symdiff::differentiate::MakeDeriv -- 277# 278# Differentiates a Tcl expression represented as a parse tree. 279# 280# Parameters: 281# tree -- Parse tree from MakeParseTreeForExpr 282# var -- Variable to differentiate with respect to 283# 284# Results: 285# Returns the parse tree of the derivative. 286 287proc math::calculus::symdiff::differentiate::MakeDeriv {tree var} { 288 return [eval [linsert $tree 1 $var]] 289} 290 291# math::calculus::symdiff::differentiate::ChainRule -- 292# 293# Applies the Chain Rule to evaluate the derivative of a unary 294# function. 295# 296# Parameters: 297# var -- Variable to differentiate with respect to. 298# derivMaker -- Command prefix for differentiating the function. 299# u -- Function argument. 300# 301# Results: 302# Returns a parse tree representing the derivative of f($u). 303# 304# ChainRule differentiates $u with respect to $var by calling MakeDeriv, 305# makes the derivative of f($u) with respect to $u by calling derivMaker 306# passing $u as a parameter, and then returns a parse tree representing 307# the product of the results. 308 309proc math::calculus::symdiff::differentiate::ChainRule {var derivMaker u} { 310 lappend derivMaker $u 311 set result [MakeProd [MakeDeriv $u $var] [eval $derivMaker]] 312} 313 314# math::calculus::symdiff::differentiate::constant -- 315# 316# Differentiate a constant. 317# 318# Parameters: 319# var -- Variable to differentiate with respect to - unused 320# constant -- Constant expression to differentiate - ignored 321# 322# Results: 323# Returns a parse tree of the derivative, which is, of course, the 324# constant zero. 325 326proc math::calculus::symdiff::differentiate::constant {var constant} { 327 return [MakeConstant 0.0] 328} 329 330# math::calculus::symdiff::differentiate::var -- 331# 332# Differentiate a variable expression. 333# 334# Parameters: 335# var - Variable with which to differentiate. 336# exprVar - Expression being differentiated, which is a single 337# variable. 338# 339# Results: 340# Returns a parse tree of the derivative. 341# 342# The derivative is the constant unity if the variables are the same 343# and the constant zero if they are different. 344 345proc math::calculus::symdiff::differentiate::var {var exprVar} { 346 if {$exprVar eq $var} { 347 return [MakeConstant 1.0] 348 } else { 349 return [MakeConstant 0.0] 350 } 351} 352 353# math::calculus::symdiff::differentiate::operator + -- 354# 355# Forms the derivative of a sum. 356# 357# Parameters: 358# var -- Variable to differentiate with respect to. 359# args -- One or two arguments giving augend and addend. If only 360# one argument is supplied, this is unary +. 361# 362# Results: 363# Returns a parse tree representing the derivative. 364# 365# Of course, the derivative of a sum is the sum of the derivatives. 366 367proc {math::calculus::symdiff::differentiate::operator +} {var args} { 368 if {[llength $args] == 1} { 369 set u [lindex $args 0] 370 set result [eval [linsert $u 1 $var]] 371 } elseif {[llength $args] == 2} { 372 foreach {u v} $args break 373 set du [eval [linsert $u 1 $var]] 374 set dv [eval [linsert $v 1 $var]] 375 set result [MakeSum $du $dv] 376 } else { 377 error "symdiff encountered a malformed parse, can't happen" 378 } 379 return $result 380} 381 382# math::calculus::symdiff::differentiate::operator - -- 383# 384# Forms the derivative of a difference. 385# 386# Parameters: 387# var -- Variable to differentiate with respect to. 388# args -- One or two arguments giving minuend and subtrahend. If only 389# one argument is supplied, this is unary -. 390# 391# Results: 392# Returns a parse tree representing the derivative. 393# 394# Of course, the derivative of a sum is the sum of the derivatives. 395 396proc {math::calculus::symdiff::differentiate::operator -} {var args} { 397 if {[llength $args] == 1} { 398 set u [lindex $args 0] 399 set du [eval [linsert $u 1 $var]] 400 set result [MakeUnaryMinus $du] 401 } elseif {[llength $args] == 2} { 402 foreach {u v} $args break 403 set du [eval [linsert $u 1 $var]] 404 set dv [eval [linsert $v 1 $var]] 405 set result [MakeDifference $du $dv] 406 } else { 407 error "symdiff encounered a malformed parse, can't happen" 408 } 409 return $result 410} 411 412# math::calculus::symdiff::differentiate::operator * -- 413# 414# Forms the derivative of a product. 415# 416# Parameters: 417# var -- Variable to differentiate with respect to. 418# u, v -- Multiplicand and multiplier. 419# 420# Results: 421# Returns a parse tree representing the derivative. 422# 423# The familiar freshman calculus product rule. 424 425proc {math::calculus::symdiff::differentiate::operator *} {var u v} { 426 set du [eval [linsert $u 1 $var]] 427 set dv [eval [linsert $v 1 $var]] 428 set result [MakeSum [MakeProd $dv $u] [MakeProd $du $v]] 429 return $result 430} 431 432# math::calculus::symdiff::differentiate::operator / -- 433# 434# Forms the derivative of a quotient. 435# 436# Parameters: 437# var -- Variable to differentiate with respect to. 438# u, v -- Dividend and divisor. 439# 440# Results: 441# Returns a parse tree representing the derivative. 442# 443# The familiar freshman calculus quotient rule. 444 445proc {math::calculus::symdiff::differentiate::operator /} {var u v} { 446 set du [eval [linsert $u 1 $var]] 447 set dv [eval [linsert $v 1 $var]] 448 set result [MakeQuotient \ 449 [MakeDifference \ 450 $du \ 451 [MakeQuotient \ 452 [MakeProd $dv $u] \ 453 $v]] \ 454 $v] 455 return $result 456} 457 458# math::calculus::symdiff::differentiate::operator acos -- 459# 460# Differentiates the 'acos' function. 461# 462# Parameters: 463# var -- Variable to differentiate with respect to. 464# u -- Argument to the acos() function. 465# 466# Results: 467# Returns a parse tree of the derivative. 468# 469# Applies the Chain Rule: D(acos(u))=-D(u)/sqrt(1 - u*u) 470# (Might it be better to factor 1-u*u into (1+u)(1-u)? Less likely to be 471# catastrophic cancellation if u is near 1?) 472 473proc {math::calculus::symdiff::differentiate::operator acos} {var u} { 474 set du [eval [linsert $u 1 $var]] 475 set result [MakeQuotient [MakeUnaryMinus $du] \ 476 [MakeFunCall sqrt \ 477 [MakeDifference [MakeConstant 1.0] \ 478 [MakeProd $u $u]]]] 479 return $result 480} 481 482# math::calculus::symdiff::differentiate::operator asin -- 483# 484# Differentiates the 'asin' function. 485# 486# Parameters: 487# var -- Variable to differentiate with respect to. 488# u -- Argument to the asin() function. 489# 490# Results: 491# Returns a parse tree of the derivative. 492# 493# Applies the Chain Rule: D(asin(u))=D(u)/sqrt(1 - u*u) 494# (Might it be better to factor 1-u*u into (1+u)(1-u)? Less likely to be 495# catastrophic cancellation if u is near 1?) 496 497proc {math::calculus::symdiff::differentiate::operator asin} {var u} { 498 set du [eval [linsert $u 1 $var]] 499 set result [MakeQuotient $du \ 500 [MakeFunCall sqrt \ 501 [MakeDifference [MakeConstant 1.0] \ 502 [MakeProd $u $u]]]] 503 return $result 504} 505 506# math::calculus::symdiff::differentiate::operator atan -- 507# 508# Differentiates the 'atan' function. 509# 510# Parameters: 511# var -- Variable to differentiate with respect to. 512# u -- Argument to the atan() function. 513# 514# Results: 515# Returns a parse tree of the derivative. 516# 517# Applies the Chain Rule: D(atan(u))=D(u)/(1 + $u*$u) 518 519proc {math::calculus::symdiff::differentiate::operator atan} {var u} { 520 set du [eval [linsert $u 1 $var]] 521 set result [MakeQuotient $du \ 522 [MakeSum [MakeConstant 1.0] \ 523 [MakeProd $u $u]]] 524} 525 526# math::calculus::symdiff::differentiate::operator atan2 -- 527# 528# Differentiates the 'atan2' function. 529# 530# Parameters: 531# var -- Variable to differentiate with respect to. 532# f, g -- Arguments to the atan() function. 533# 534# Results: 535# Returns a parse tree of the derivative. 536# 537# Applies the Chain and Quotient Rules: 538# D(atan2(f, g)) = (D(f)*g - D(g)*f)/(f*f + g*g) 539 540proc {math::calculus::symdiff::differentiate::operator atan2} {var f g} { 541 set df [eval [linsert $f 1 $var]] 542 set dg [eval [linsert $g 1 $var]] 543 return [MakeQuotient \ 544 [MakeDifference \ 545 [MakeProd $df $g] \ 546 [MakeProd $f $dg]] \ 547 [MakeSum \ 548 [MakeProd $f $f] \ 549 [MakeProd $g $g]]] 550} 551 552# math::calculus::symdiff::differentiate::operator cos -- 553# 554# Differentiates the 'cos' function. 555# 556# Parameters: 557# var -- Variable to differentiate with respect to. 558# u -- Argument to the cos() function. 559# 560# Results: 561# Returns a parse tree of the derivative. 562# 563# Applies the Chain Rule: D(cos(u))=-sin(u)*D(u) 564 565proc {math::calculus::symdiff::differentiate::operator cos} {var u} { 566 return [ChainRule $var MakeMinusSin $u] 567} 568proc math::calculus::symdiff::differentiate::MakeMinusSin {operand} { 569 return [MakeUnaryMinus [MakeFunCall sin $operand]] 570} 571 572# math::calculus::symdiff::differentiate::operator cosh -- 573# 574# Differentiates the 'cosh' function. 575# 576# Parameters: 577# var -- Variable to differentiate with respect to. 578# u -- Argument to the cosh() function. 579# 580# Results: 581# Returns a parse tree of the derivative. 582# 583# Applies the Chain Rule: D(cosh(u))=sinh(u)*D(u) 584 585proc {math::calculus::symdiff::differentiate::operator cosh} {var u} { 586 set result [ChainRule $var [list MakeFunCall sinh] $u] 587 return $result 588} 589 590# math::calculus::symdiff::differentiate::operator exp -- 591# 592# Differentiate the exponential function 593# 594# Parameters: 595# var -- Variable to differentiate with respect to. 596# u -- Argument of the exponential function. 597# 598# Results: 599# Returns a parse tree of the derivative. 600# 601# Uses the Chain Rule D(exp(u)) = exp(u)*D(u). 602 603proc {math::calculus::symdiff::differentiate::operator exp} {var u} { 604 set result [ChainRule $var [list MakeFunCall exp] $u] 605 return $result 606} 607 608# math::calculus::symdiff::differentiate::operator hypot -- 609# 610# Differentiate the 'hypot' function 611# 612# Parameters: 613# var - Variable to differentiate with respect to. 614# f, g - Arguments to the 'hypot' function 615# 616# Results: 617# Returns a parse tree of the derivative 618# 619# Uses a number of algebraic simplifications to arrive at: 620# D(hypot(f,g)) = (f*D(f)+g*D(g))/hypot(f,g) 621 622proc {math::calculus::symdiff::differentiate::operator hypot} {var f g} { 623 set df [eval [linsert $f 1 $var]] 624 set dg [eval [linsert $g 1 $var]] 625 return [MakeQuotient \ 626 [MakeSum \ 627 [MakeProd $df $f] \ 628 [MakeProd $dg $g]] \ 629 [MakeFunCall hypot $f $g]] 630} 631 632# math::calculus::symdiff::differentiate::operator log -- 633# 634# Differentiates a logarithm. 635# 636# Parameters: 637# var -- Variable to differentiate with respect to. 638# u -- Argument to the log() function. 639# 640# Results: 641# Returns a parse tree of the derivative. 642# 643# D(log(u))==D(u)/u 644 645proc {math::calculus::symdiff::differentiate::operator log} {var u} { 646 set du [eval [linsert $u 1 $var]] 647 set result [MakeQuotient $du $u] 648 return $result 649} 650 651# math::calculus::symdiff::differentiate::operator log10 -- 652# 653# Differentiates a common logarithm. 654# 655# Parameters: 656# var -- Variable to differentiate with respect to. 657# u -- Argument to the log10() function. 658# 659# Results: 660# Returns a parse tree of the derivative. 661# 662# D(log(u))==D(u)/(u * log(10)) 663 664proc {math::calculus::symdiff::differentiate::operator log10} {var u} { 665 set du [eval [linsert $u 1 $var]] 666 set result [MakeQuotient $du \ 667 [MakeProd [MakeConstant [expr log(10.)]] $u]] 668 return $result 669} 670 671# math::calculus::symdiff::differentiate::operator ** -- 672# 673# Differentiate an exponential. 674# 675# Parameters: 676# var -- Variable to differentiate with respect to 677# f, g -- Base and exponent 678# 679# Results: 680# Returns the parse tree of the derivative. 681# 682# Handles the special case where g is constant as 683# D(f**g) == g*f**(g-1)*D(f) 684# Otherwise, uses the general power formula 685# D(f**g) == (f**g) * (((D(f)*g)/f) + (D(g)*log(f))) 686 687proc {math::calculus::symdiff::differentiate::operator **} {var f g} { 688 set df [eval [linsert $f 1 $var]] 689 if {[IsConstant $g]} { 690 set gm1 [MakeConstant [expr {[ConstantValue $g] - 1}]] 691 set result [MakeProd $df [MakeProd $g [MakePower $f $gm1]]] 692 693 } else { 694 set dg [eval [linsert $g 1 $var]] 695 set result [MakeProd [MakePower $f $g] \ 696 [MakeSum \ 697 [MakeQuotient [MakeProd $df $g] $f] \ 698 [MakeProd $dg [MakeFunCall log $f]]]] 699 } 700 return $result 701} 702interp alias {} {math::calculus::symdiff::differentiate::operator pow} \ 703 {} {math::calculus::symdiff::differentiate::operator **} 704 705# math::calculus::symdiff::differentiate::operator sin -- 706# 707# Differentiates the 'sin' function. 708# 709# Parameters: 710# var -- Variable to differentiate with respect to. 711# u -- Argument to the sin() function. 712# 713# Results: 714# Returns a parse tree of the derivative. 715# 716# Applies the Chain Rule: D(sin(u))=cos(u)*D(u) 717 718proc {math::calculus::symdiff::differentiate::operator sin} {var u} { 719 set result [ChainRule $var [list MakeFunCall cos] $u] 720 return $result 721} 722 723# math::calculus::symdiff::differentiate::operator sinh -- 724# 725# Differentiates the 'sinh' function. 726# 727# Parameters: 728# var -- Variable to differentiate with respect to. 729# u -- Argument to the sinh() function. 730# 731# Results: 732# Returns a parse tree of the derivative. 733# 734# Applies the Chain Rule: D(sin(u))=cosh(u)*D(u) 735 736proc {math::calculus::symdiff::differentiate::operator sinh} {var u} { 737 set result [ChainRule $var [list MakeFunCall cosh] $u] 738 return $result 739} 740 741# math::calculus::symdiff::differentiate::operator sqrt -- 742# 743# Differentiate the 'sqrt' function. 744# 745# Parameters: 746# var -- Variable to differentiate with respect to 747# u -- Parameter of 'sqrt' as a parse tree. 748# 749# Results: 750# Returns a parse tree representing the derivative. 751# 752# D(sqrt(u))==D(u)/(2*sqrt(u)) 753 754proc {math::calculus::symdiff::differentiate::operator sqrt} {var u} { 755 set du [eval [linsert $u 1 $var]] 756 set result [MakeQuotient $du [MakeProd [MakeConstant 2.0] \ 757 [MakeFunCall sqrt $u]]] 758 return $result 759} 760 761# math::calculus::symdiff::differentiate::operator tan -- 762# 763# Differentiates the 'tan' function. 764# 765# Parameters: 766# var -- Variable to differentiate with respect to. 767# u -- Argument to the tan() function. 768# 769# Results: 770# Returns a parse tree of the derivative. 771# 772# Applies the Chain Rule: D(tan(u))=D(u)/(cos(u)*cos(u)) 773 774proc {math::calculus::symdiff::differentiate::operator tan} {var u} { 775 set du [eval [linsert $u 1 $var]] 776 set cosu [MakeFunCall cos $u] 777 return [MakeQuotient $du [MakeProd $cosu $cosu]] 778} 779 780# math::calculus::symdiff::differentiate::operator tanh -- 781# 782# Differentiates the 'tanh' function. 783# 784# Parameters: 785# var -- Variable to differentiate with respect to. 786# u -- Argument to the tanh() function. 787# 788# Results: 789# Returns a parse tree of the derivative. 790# 791# Applies the Chain Rule: D(tanh(u))=D(u)/(cosh(u)*cosh(u)) 792 793proc {math::calculus::symdiff::differentiate::operator tanh} {var u} { 794 set du [eval [linsert $u 1 $var]] 795 set coshu [MakeFunCall cosh $u] 796 return [MakeQuotient $du [MakeProd $coshu $coshu]] 797} 798 799# math::calculus::symdiff::MakeFunCall -- 800# 801# Makes a parse tree for a function call 802# 803# Parameters: 804# fun -- Name of the function to call 805# args -- Arguments to the function, expressed as parse trees 806# 807# Results: 808# Returns a parse tree for the result of calling the function. 809# 810# Performs the peephole optimization of replacing a function with 811# constant parameters with its value. 812 813proc math::calculus::symdiff::MakeFunCall {fun args} { 814 set constant 1 815 set exp $fun 816 append exp \( 817 set sep "" 818 foreach a $args { 819 if {[IsConstant $a]} { 820 append exp $sep [ConstantValue $a] 821 set sep "," 822 } else { 823 set constant 0 824 break 825 } 826 } 827 if {$constant} { 828 append exp \) 829 return [MakeConstant [expr $exp]] 830 } 831 set result [MakeOperator $fun] 832 foreach arg $args { 833 lappend result $arg 834 } 835 return $result 836} 837 838# math::calculus::symdiff::MakeSum -- 839# 840# Makes the parse tree for a sum. 841# 842# Parameters: 843# left, right -- Parse trees for augend and addend 844# 845# Results: 846# Returns the parse tree for the sum. 847# 848# Performs the following peephole optimizations: 849# (1) a + (-b) = a - b 850# (2) (-a) + b = b - a 851# (3) 0 + a = a 852# (4) a + 0 = a 853# (5) The sum of two constants may be reduced to a constant 854 855proc math::calculus::symdiff::MakeSum {left right} { 856 if {[IsUnaryMinus $right]} { 857 return [MakeDifference $left [UnaryMinusArg $right]] 858 } 859 if {[IsUnaryMinus $left]} { 860 return [MakeDifference $right [UnaryMinusArg $left]] 861 } 862 if {[IsConstant $left]} { 863 set v [ConstantValue $left] 864 if {$v == 0} { 865 return $right 866 } elseif {[IsConstant $right]} { 867 return [MakeConstant [expr {[ConstantValue $left] 868 + [ConstantValue $right]}]] 869 } 870 } elseif {[IsConstant $right]} { 871 set v [ConstantValue $right] 872 if {$v == 0} { 873 return $left 874 } 875 } 876 set result [MakeOperator +] 877 lappend result $left $right 878 return $result 879} 880 881# math::calculus::symdiff::MakeDifference -- 882# 883# Makes the parse tree for a difference 884# 885# Parameters: 886# left, right -- Minuend and subtrahend, expressed as parse trees 887# 888# Results: 889# Returns a parse tree expressing the difference 890# 891# Performs the following peephole optimizations: 892# (1) a - (-b) = a + b 893# (2) -a - b = -(a + b) 894# (3) 0 - b = -b 895# (4) a - 0 = a 896# (5) The difference of any two constants can be reduced to a constant. 897 898proc math::calculus::symdiff::MakeDifference {left right} { 899 if {[IsUnaryMinus $right]} { 900 return [MakeSum $left [UnaryMinusArg $right]] 901 } 902 if {[IsUnaryMinus $left]} { 903 return [MakeUnaryMinus [MakeSum [UnaryMinusArg $left] $right]] 904 } 905 if {[IsConstant $left]} { 906 set v [ConstantValue $left] 907 if {$v == 0} { 908 return [MakeUnaryMinus $right] 909 } elseif {[IsConstant $right]} { 910 return [MakeConstant [expr {[ConstantValue $left] 911 - [ConstantValue $right]}]] 912 } 913 } elseif {[IsConstant $right]} { 914 set v [ConstantValue $right] 915 if {$v == 0} { 916 return $left 917 } 918 } 919 set result [MakeOperator -] 920 lappend result $left $right 921 return $result 922} 923 924# math::calculus::symdiff::MakeProd -- 925# 926# Constructs the parse tree for a product, left*right. 927# 928# Parameters: 929# left, right - Multiplicand and multiplier 930# 931# Results: 932# Returns the parse tree for the result. 933# 934# Performs the following peephole optimizations. 935# (1) If either operand is a unary minus, it is hoisted out of the 936# expression. 937# (2) If either operand is the constant 0, the result is the constant 0 938# (3) If either operand is the constant 1, the result is the other operand. 939# (4) If either operand is the constant -1, the result is unary minus 940# applied to the other operand 941# (5) If both operands are constant, the result is a constant containing 942# their product. 943 944proc math::calculus::symdiff::MakeProd {left right} { 945 if {[IsUnaryMinus $left]} { 946 return [MakeUnaryMinus [MakeProd [UnaryMinusArg $left] $right]] 947 } 948 if {[IsUnaryMinus $right]} { 949 return [MakeUnaryMinus [MakeProd $left [UnaryMinusArg $right]]] 950 } 951 if {[IsConstant $left]} { 952 set v [ConstantValue $left] 953 if {$v == 0} { 954 return [MakeConstant 0.0] 955 } elseif {$v == 1} { 956 return $right 957 } elseif {$v == -1} { 958 return [MakeUnaryMinus $right] 959 } elseif {[IsConstant $right]} { 960 return [MakeConstant [expr {[ConstantValue $left] 961 * [ConstantValue $right]}]] 962 } 963 } elseif {[IsConstant $right]} { 964 set v [ConstantValue $right] 965 if {$v == 0} { 966 return [MakeConstant 0.0] 967 } elseif {$v == 1} { 968 return $left 969 } elseif {$v == -1} { 970 return [MakeUnaryMinus $left] 971 } 972 } 973 set result [MakeOperator *] 974 lappend result $left $right 975 return $result 976} 977 978# math::calculus::symdiff::MakeQuotient -- 979# 980# Makes a parse tree for a quotient, n/d 981# 982# Parameters: 983# n, d - Parse trees for numerator and denominator 984# 985# Results: 986# Returns the parse tree for the quotient. 987# 988# Performs peephole optimizations: 989# (1) If either operand is a unary minus, it is hoisted out. 990# (2) If the numerator is the constant 0, the result is the constant 0. 991# (3) If the demominator is the constant 1, the result is the numerator 992# (4) If the denominator is the constant -1, the result is the unary 993# negation of the numerator. 994# (5) If both numerator and denominator are constant, the result is 995# a constant representing their quotient. 996 997proc math::calculus::symdiff::MakeQuotient {n d} { 998 if {[IsUnaryMinus $n]} { 999 return [MakeUnaryMinus [MakeQuotient [UnaryMinusArg $n] $d]] 1000 } 1001 if {[IsUnaryMinus $d]} { 1002 return [MakeUnaryMinus [MakeQuotient $n [UnaryMinusArg $d]]] 1003 } 1004 if {[IsConstant $n]} { 1005 set v [ConstantValue $n] 1006 if {$v == 0} { 1007 return [MakeConstant 0.0] 1008 } elseif {[IsConstant $d]} { 1009 return [MakeConstant [expr {[ConstantValue $n] 1010 * [ConstantValue $d]}]] 1011 } 1012 } elseif {[IsConstant $d]} { 1013 set v [ConstantValue $d] 1014 if {$v == 0} { 1015 return -code error "requested expression will result in division by zero at run time" 1016 } elseif {$v == 1} { 1017 return $n 1018 } elseif {$v == -1} { 1019 return [MakeUnaryMinus $n] 1020 } 1021 } 1022 set result [MakeOperator /] 1023 lappend result $n $d 1024 return $result 1025} 1026 1027# math::calculus::symdiff::MakePower -- 1028# 1029# Make a parse tree for an exponentiation operation 1030# 1031# Parameters: 1032# a -- Base, expressed as a parse tree 1033# b -- Exponent, expressed as a parse tree 1034# 1035# Results: 1036# Returns a parse tree for the expression 1037# 1038# Performs peephole optimizations: 1039# (1) The constant zero raised to any non-zero power is 0 1040# (2) The constant 1 raised to any power is 1 1041# (3) Any non-zero quantity raised to the zero power is 1 1042# (4) Any non-zero quantity raised to the first power is the base itself. 1043# (5) MakeFunCall will optimize any other case of a constant raised 1044# to a constant power. 1045 1046proc math::calculus::symdiff::MakePower {a b} { 1047 if {[IsConstant $a]} { 1048 if {[ConstantValue $a] == 0} { 1049 if {[IsConstant $b] && [ConstantValue $b] == 0} { 1050 error "requested expression will result in zero to zero power at run time" 1051 } 1052 return [MakeConstant 0.0] 1053 } elseif {[ConstantValue $a] == 1} { 1054 return [MakeConstant 1.0] 1055 } 1056 } 1057 if {[IsConstant $b]} { 1058 if {[ConstantValue $b] == 0} { 1059 return [MakeConstant 1.0] 1060 } elseif {[ConstantValue $b] == 1} { 1061 return $a 1062 } 1063 } 1064 return [MakeFunCall pow $a $b] 1065} 1066 1067# math::calculus::symdiff::MakeUnaryMinus -- 1068# 1069# Makes the parse tree for a unary negation. 1070# 1071# Parameters: 1072# operand -- Parse tree for the operand 1073# 1074# Results: 1075# Returns the parse tree for the expression 1076# 1077# Performs the following peephole optimizations: 1078# (1) -(-$a) = $a 1079# (2) The unary negation of a constant is another constant 1080 1081proc math::calculus::symdiff::MakeUnaryMinus {operand} { 1082 if {[IsUnaryMinus $operand]} { 1083 return [UnaryMinusArg $operand] 1084 } 1085 if {[IsConstant $operand]} { 1086 return [MakeConstant [expr {-[ConstantValue $operand]}]] 1087 } else { 1088 return [list [list operator -] $operand] 1089 } 1090} 1091 1092# math::calculus::symdiff::IsUnaryMinus -- 1093# 1094# Determines whether a parse tree represents a unary negation 1095# 1096# Parameters: 1097# x - Parse tree to examine 1098# 1099# Results: 1100# Returns 1 if the parse tree represents a unary minus, 0 otherwise 1101 1102proc math::calculus::symdiff::IsUnaryMinus {x} { 1103 return [expr {[llength $x] == 2 1104 && [lindex $x 0] eq [list operator -]}] 1105} 1106 1107# math::calculus::symdiff::UnaryMinusArg -- 1108# 1109# Extracts the argument from a unary negation. 1110# 1111# Parameters: 1112# x - Parse tree to examine, known to represent a unary negation 1113# 1114# Results: 1115# Returns a parse tree representing the operand. 1116 1117proc math::calculus::symdiff::UnaryMinusArg {x} { 1118 return [lindex $x 1] 1119} 1120 1121# math::calculus::symdiff::MakeOperator -- 1122# 1123# Makes a partial parse tree for an operator 1124# 1125# Parameters: 1126# op -- Name of the operator 1127# 1128# Results: 1129# Returns the resulting parse tree. 1130# 1131# The caller may use [lappend] to place any needed operands 1132 1133proc math::calculus::symdiff::MakeOperator {op} { 1134 if {$op eq {?}} { 1135 return -code error "symdiff can't differentiate the ternary ?: operator" 1136 } elseif {[namespace which [list differentiate::operator $op]] ne {}} { 1137 return [list [list operator $op]] 1138 } elseif {[string is alnum $op] && ($op ni {eq ne in ni})} { 1139 return -code error "symdiff can't differentiate the \"$op\" function" 1140 } else { 1141 return -code error "symdiff can't differentiate the \"$op\" operator" 1142 } 1143} 1144 1145# math::calculus::symdiff::MakeVariable -- 1146# 1147# Makes a partial parse tree for a single variable 1148# 1149# Parameters: 1150# name -- Name of the variable 1151# 1152# Results: 1153# Returns a partial parse tree giving the variable 1154 1155proc math::calculus::symdiff::MakeVariable {name} { 1156 return [list var $name] 1157} 1158 1159# math::calculus::symdiff::MakeConstant -- 1160# 1161# Make the parse tree for a constant. 1162# 1163# Parameters: 1164# value -- The constant's value 1165# 1166# Results: 1167# Returns a parse tree. 1168 1169proc math::calculus::symdiff::MakeConstant {value} { 1170 return [list constant $value] 1171} 1172 1173# math::calculus::symdiff::IsConstant -- 1174# 1175# Test if an expression represented by a parse tree is a constant. 1176# 1177# Parameters: 1178# Item - Parse tree to test 1179# 1180# Results: 1181# Returns 1 for a constant, 0 for anything else 1182 1183proc math::calculus::symdiff::IsConstant {item} { 1184 return [expr {[lindex $item 0] eq {constant}}] 1185} 1186 1187# math::calculus::symdiff::ConstantValue -- 1188# 1189# Recovers a constant value from the parse tree representing a constant 1190# expression. 1191# 1192# Parameters: 1193# item -- Parse tree known to be a constant. 1194# 1195# Results: 1196# Returns the constant value. 1197 1198proc math::calculus::symdiff::ConstantValue {item} { 1199 return [lindex $item 1] 1200} 1201 1202# Define the parse tree fabrication routines in the 'differentiate' 1203# namespace as well as the 'symdiff' namespace, without exporting them 1204# from the package. 1205 1206interp alias {} math::calculus::symdiff::differentiate::IsConstant \ 1207 {} math::calculus::symdiff::IsConstant 1208interp alias {} math::calculus::symdiff::differentiate::ConstantValue \ 1209 {} math::calculus::symdiff::ConstantValue 1210interp alias {} math::calculus::symdiff::differentiate::MakeConstant \ 1211 {} math::calculus::symdiff::MakeConstant 1212interp alias {} math::calculus::symdiff::differentiate::MakeDifference \ 1213 {} math::calculus::symdiff::MakeDifference 1214interp alias {} math::calculus::symdiff::differentiate::MakeFunCall \ 1215 {} math::calculus::symdiff::MakeFunCall 1216interp alias {} math::calculus::symdiff::differentiate::MakePower \ 1217 {} math::calculus::symdiff::MakePower 1218interp alias {} math::calculus::symdiff::differentiate::MakeProd \ 1219 {} math::calculus::symdiff::MakeProd 1220interp alias {} math::calculus::symdiff::differentiate::MakeQuotient \ 1221 {} math::calculus::symdiff::MakeQuotient 1222interp alias {} math::calculus::symdiff::differentiate::MakeSum \ 1223 {} math::calculus::symdiff::MakeSum 1224interp alias {} math::calculus::symdiff::differentiate::MakeUnaryMinus \ 1225 {} math::calculus::symdiff::MakeUnaryMinus 1226interp alias {} math::calculus::symdiff::differentiate::MakeVariable \ 1227 {} math::calculus::symdiff::MakeVariable 1228interp alias {} math::calculus::symdiff::differentiate::ExtractExpression \ 1229 {} math::calculus::symdiff::ExtractExpression 1230