Merge commit '7c621da28f07e449ad90c387344f07a453927569' as 'tcllib'

1 files changed, 364 insertions, 0 deletions

diff --git a/tcllib/modules/math/rational_funcs.tcl b/tcllib/modules/math/rational_funcs.tcl
new file mode 100755
index 0000000..2dae397
--- /dev/null
+++ b/tcllib/modules/math/rational_funcs.tcl
@@ -0,0 +1,364 @@
+# rational_funcs.tcl --
+#    Implement procedures to deal with rational functions
+#
+
+package require math::polynomials
+
+namespace eval ::math::rationalfunctions {
+    variable count 0  ;# Count the number of specific commands
+    namespace eval v {}
+
+    namespace export rationalFunction ratioCmd evalRatio \
+                     coeffsNumerator coeffsDenominator \
+                     derivRatio  \
+                     addRatio    subRatio multRatio \
+                     divRatio
+
+    namespace import ::math::polynomials::*
+}
+
+
+# rationalFunction --
+#    Return a rational function definition
+#
+# Arguments:
+#    num          The coefficients of the numerator
+#    den          The coefficients of the denominator
+# Result:
+#    Rational function definition
+#
+proc ::math::rationalfunctions::rationalFunction {num den} {
+
+    foreach coeffs [list $num $den] {
+        foreach coeff $coeffs {
+            if { ! [string is double -strict $coeff] } {
+                return -code error "Coefficients must be real numbers"
+            }
+        }
+    }
+
+    #
+    # The leading coefficient must be non-zero
+    #
+    return [list RATIONAL_FUNCTION [polynomial $num] [polynomial $den]]
+}
+
+# ratioCmd --
+#    Return a procedure that implements a rational function evaluation
+#
+# Arguments:
+#    num          The coefficients of the numerator
+#    den          The coefficients of the denominator
+# Result:
+#    New procedure
+#
+proc ::math::rationalfunctions::ratioCmd {num {den {}}} {
+    variable count
+
+    if { [llength $den] == 0 } {
+        if { [lindex $num 0] == "RATIONAL_FUNCTION" } {
+            set den [lindex $num 2]
+            set num [lindex $num 1]
+        }
+    }
+
+    set degree1 [expr {[llength $num]-1}]
+    set degree2 [expr {[llength $num]-1}]
+    set body "expr \{([join $num +\$x*(][string repeat ) $degree1])/\
+(double([join $den +\$x*(][string repeat ) $degree2])\}"
+
+    incr count
+    set name "::math::rationalfunctions::v::RATIO$count"
+    proc $name {x} $body
+    return $name
+}
+
+# evalRatio --
+#    Evaluate a rational function at a given coordinate
+#
+# Arguments:
+#    ratio        Rational function definition
+#    x            Coordinate
+# Result:
+#    Value at x
+#
+proc ::math::rationalfunctions::evalRatio {ratio x} {
+    if { [lindex $ratio 0] != "RATIONAL_FUNCTION" } {
+        return -code error "Not a rational function"
+    }
+    if { ! [string is double $x] } {
+        return -code error "Coordinate must be a real number"
+    }
+
+    set num 0.0
+    foreach c [lindex [lindex $ratio 1] 1] {
+        set num [expr {$num*$x+$c}]
+    }
+
+    set den 0.0
+    foreach c [lindex [lindex $ratio 2] 1] {
+        set den [expr {$den*$x+$c}]
+    }
+    return [expr {$num/double($den)}]
+}
+
+# coeffsNumerator --
+#    Return the coefficients of the numerator
+#
+# Arguments:
+#    ratio        Rational function definition
+# Result:
+#    The coefficients in ascending order
+#
+proc ::math::rationalfunctions::coeffsNumerator {ratio} {
+    if { [lindex $ratio 0] != "RATIONAL_FUNCTION" } {
+        return -code error "Not a rational function"
+    }
+    set polyn [lindex $ratio 1]
+    return [allCoeffsPolyn $polyn]
+}
+
+# coeffsDenominator --
+#    Return the coefficients of the denominator
+#
+# Arguments:
+#    ratio        Rational function definition
+# Result:
+#    The coefficients in ascending order
+#
+proc ::math::rationalfunctions::coeffsDenominator {ratio} {
+    if { [lindex $ratio 0] != "RATIONAL_FUNCTION" } {
+        return -code error "Not a rational function"
+    }
+    set polyn [lindex $ratio 2]
+    return [allCoeffsPolyn $polyn]
+}
+
+# derivRatio --
+#    Return the derivative of the rational function
+#
+# Arguments:
+#    polyn        Polynomial definition
+# Result:
+#    The new polynomial
+#
+proc ::math::rationalfunctions::derivRatio {ratio} {
+    if { [lindex $ratio 0] != "RATIONAL_FUNCTION" } {
+        return -code error "Not a rational function"
+    }
+    set num_polyn [lindex $ratio 1]
+    set den_polyn [lindex $ratio 2]
+    set num_deriv [derivPolyn $num_polyn]
+    set den_deriv [derivPolyn $den_polyn]
+    set num       [subPolyn [multPolyn $num_deriv $den_polyn] \
+                            [multPolyn $den_deriv $num_polyn] ]
+    set den       [multPolyn $den_polyn $den_polyn]
+
+    return [list RATIONAL_FUNCTION $num $den]
+}
+
+# addRatio --
+#    Add two rational functions and return the result
+#
+# Arguments:
+#    ratio1       First rational function or a scalar
+#    ratio2       Second rational function or a scalar
+# Result:
+#    The sum of the two functions
+# Note:
+#    TODO: Check for the same denominator
+#
+proc ::math::rationalfunctions::addRatio {ratio1 ratio2} {
+    if { [llength $ratio1] == 1 && [string is double -strict $ratio1] } {
+        set polyn1 [rationalFunction $ratio1 1.0]
+    }
+    if { [llength $ratio2] == 1 && [string is double -strict $ratio2] } {
+        set ratio2 [rationalFunction $ratio1 1.0]
+    }
+    if { [lindex $ratio1 0] != "RATIONAL_FUNCTION" ||
+         [lindex $ratio2 0] != "RATIONAL_FUNCTION" } {
+        return -code error "Both arguments must be rational functions or a real number"
+    }
+
+    set num1    [lindex $ratio1 1]
+    set den1    [lindex $ratio1 2]
+    set num2    [lindex $ratio2 1]
+    set den2    [lindex $ratio2 2]
+
+    set newnum  [addPolyn [multPolyn $num1 $den2] \
+                          [multPolyn $num2 $den1] ]
+
+    set newden  [multPolyn $den1 $den2]
+
+    return [list RATIONAL_FUNCTION $newnum $newden]
+}
+
+# subRatio --
+#    Subtract two rational functions and return the result
+#
+# Arguments:
+#    ratio1       First rational function or a scalar
+#    ratio2       Second rational function or a scalar
+# Result:
+#    The difference of the two functions
+# Note:
+#    TODO: Check for the same denominator
+#
+proc ::math::rationalfunctions::subRatio {ratio1 ratio2} {
+    if { [llength $ratio1] == 1 && [string is double -strict $ratio1] } {
+        set polyn1 [rationalFunction $ratio1 1.0]
+    }
+    if { [llength $ratio2] == 1 && [string is double -strict $ratio2] } {
+        set ratio2 [rationalFunction $ratio1 1.0]
+    }
+    if { [lindex $ratio1 0] != "RATIONAL_FUNCTION" ||
+         [lindex $ratio2 0] != "RATIONAL_FUNCTION" } {
+        return -code error "Both arguments must be rational functions or a real number"
+    }
+
+    set num1    [lindex $ratio1 1]
+    set den1    [lindex $ratio1 2]
+    set num2    [lindex $ratio2 1]
+    set den2    [lindex $ratio2 2]
+
+    set newnum  [subPolyn [multPolyn $num1 $den2] \
+                          [multPolyn $num2 $den1] ]
+
+    set newden  [multPolyn $den1 $den2]
+
+    return [list RATIONAL_FUNCTION $newnum $newden]
+}
+
+# multRatio --
+#    Multiply two rational functions and return the result
+#
+# Arguments:
+#    ratio1       First rational function or a scalar
+#    ratio2       Second rational function or a scalar
+# Result:
+#    The product of the two functions
+# Note:
+#    TODO: Check for the same denominator
+#
+proc ::math::rationalfunctions::multRatio {ratio1 ratio2} {
+    if { [llength $ratio1] == 1 && [string is double -strict $ratio1] } {
+        set polyn1 [rationalFunction $ratio1 1.0]
+    }
+    if { [llength $ratio2] == 1 && [string is double -strict $ratio2] } {
+        set ratio2 [rationalFunction $ratio1 1.0]
+    }
+    if { [lindex $ratio1 0] != "RATIONAL_FUNCTION" ||
+         [lindex $ratio2 0] != "RATIONAL_FUNCTION" } {
+        return -code error "Both arguments must be rational functions or a real number"
+    }
+
+    set num1    [lindex $ratio1 1]
+    set den1    [lindex $ratio1 2]
+    set num2    [lindex $ratio2 1]
+    set den2    [lindex $ratio2 2]
+
+    set newnum  [multPolyn $num1 $num2]
+    set newden  [multPolyn $den1 $den2]
+
+    return [list RATIONAL_FUNCTION $newnum $newden]
+}
+
+# divRatio --
+#    Divide two rational functions and return the result
+#
+# Arguments:
+#    ratio1       First rational function or a scalar
+#    ratio2       Second rational function or a scalar
+# Result:
+#    The quotient of the two functions
+# Note:
+#    TODO: Check for the same denominator
+#
+proc ::math::rationalfunctions::divRatio {ratio1 ratio2} {
+    if { [llength $ratio1] == 1 && [string is double -strict $ratio1] } {
+        set polyn1 [rationalFunction $ratio1 1.0]
+    }
+    if { [llength $ratio2] == 1 && [string is double -strict $ratio2] } {
+        set ratio2 [rationalFunction $ratio1 1.0]
+    }
+    if { [lindex $ratio1 0] != "RATIONAL_FUNCTION" ||
+         [lindex $ratio2 0] != "RATIONAL_FUNCTION" } {
+        return -code error "Both arguments must be rational functions or a real number"
+    }
+
+    set num1    [lindex $ratio1 1]
+    set den1    [lindex $ratio1 2]
+    set num2    [lindex $ratio2 1]
+    set den2    [lindex $ratio2 2]
+
+    set newnum  [multPolyn $num1 $den2]
+    set newden  [multPolyn $num2 $den1]
+
+    return [list RATIONAL_FUNCTION $newnum $newden]
+}
+
+#
+# Announce our presence
+#
+package provide math::rationalfunctions 1.0.1
+
+# some tests --
+#
+if { 0 } {
+set prec $::tcl_precision
+if {![package vsatisfies [package provide Tcl] 8.5]} {
+    set ::tcl_precision 17
+} else {
+    set ::tcl_precision 0
+}
+
+
+set f1    [::math::rationalfunctions::rationalFunction {1 2 3} {1 4}]
+set f2    [::math::rationalfunctions::rationalFunction {1 2 3 0} {1 4}]
+set f3    [::math::rationalfunctions::rationalFunction {0 0 0 0} {1}]
+set f4    [::math::rationalfunctions::rationalFunction {5 7} {1}]
+set cmdf1 [::math::rationalfunctions::ratioCmd {1 2 3} {1 4}]
+
+foreach x {0 1 2 3 4 5} {
+    puts "[::math::rationalfunctions::evalRatio $f1 $x] -- \
+[expr {(1.0+2.0*$x+3.0*$x*$x)/double(1.0+4.0*$x)}] -- \
+[$cmdf1 $x] -- [::math::rationalfunctions::evalRatio $f3 $x]"
+}
+
+puts "All coefficients = [::math::rationalfunctions::coeffsNumerator $f2]"
+puts "                   [::math::rationalfunctions::coeffsDenominator $f2]"
+
+puts "Derivative = [::math::rationalfunctions::derivRatio $f1]"
+
+puts "Add:       [::math::rationalfunctions::addRatio $f1 $f4]"
+puts "Add:       [::math::rationalfunctions::addRatio $f4 $f1]"
+puts "Subtract:  [::math::rationalfunctions::subRatio $f1 $f4]"
+puts "Multiply:  [::math::rationalfunctions::multRatio $f1 $f4]"
+
+set f1    [::math::rationalfunctions::rationalFunction {1 2 3} 1]
+set f2    [::math::rationalfunctions::rationalFunction {0 1} 1]
+
+puts "Divide:    [::math::rationalfunctions::divRatio $f1 $f2]"
+
+set f1    [::math::rationalfunctions::rationalFunction {1 2 3} 1]
+set f2    [::math::rationalfunctions::rationalFunction {1 1} {1 2}]
+
+puts "Divide:    [::math::rationalfunctions::divRatio $f1 $f2]"
+
+set f1 [::math::rationalfunctions::rationalFunction {1 2 3} 1]
+set f2 [::math::rationalfunctions::rationalFunction {0 1} {0 0 1}]
+set f3 [::math::rationalfunctions::divRatio $f2 $f1]
+set coeffs [::math::rationalfunctions::coeffsNumerator $f3]
+puts "Coefficients: $coeffs"
+set f3 [::math::rationalfunctions::divRatio $f1 $f2]
+set coeffs [::math::rationalfunctions::coeffsNumerator $f3]
+puts "Coefficients: $coeffs"
+set f1 [::math::rationalfunctions::rationalFunction {1 2 3} {1 2}]
+set f2 [::math::rationalfunctions::rationalFunction {0} {1}]
+set f3 [::math::rationalfunctions::divRatio $f2 $f1]
+set coeffs [::math::rationalfunctions::coeffsNumerator $f3]
+puts "Coefficients: $coeffs"
+puts "Eval null function: [::math::rationalfunctions::evalRatio $f2 1]"
+
+set ::tcl_precision $prec
+}