rational_funcs.tcl

Go to the documentation of this file.

/*  rational_funcs.tcl --*/
/*     Implement procedures to deal with rational functions*/
/* */

package require math::polynomials

namespace ::math::rationalfunctions {
    variable count 0  ;/*  Count the number of specific commands*/
    namespace 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*/
/* */
ret  ::math::rationalfunctions::rationalFunction (type num , type 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*/
/* */
ret  ::math::rationalfunctions::ratioCmd (type num , optional den ={)} {
    variable count

    if { [llength $den] == 0 } {
        if { [lindex $num 0] == "RATIONAL_FUNCTION" } {
             den =  [lindex $num 2]
             num =  [lindex $num 1]
        }
    }

     degree1 =  [expr {[llength $num]-1}]
     degree2 =  [expr {[llength $num]-1}]
     body =  "expr \{([join $num +\$x*(][string repeat ) $degree1])/\
(double([join $den +\$x*(][string repeat ) $degree2])\}"

    incr count
     name =  "::math::rationalfunctions::v::RATIO$count"
    ret  $name (type 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"
    }

     num =  0.0
    foreach c [lindex [lindex $ratio 1] 1] {
         num =  [expr {$num*$x+$c}]
    }

     den =  0.0
    foreach c [lindex [lindex $ratio 2] 1] {
         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*/
/* */
ret  ::math::rationalfunctions::coeffsNumerator (type 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*/
/* */
ret  ::math::rationalfunctions::coeffsDenominator (type 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*/
/* */
ret  ::math::rationalfunctions::derivRatio (type 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*/
/* */
ret  ::math::rationalfunctions::addRatio (type ratio1 , type 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*/
/* */
ret  ::math::rationalfunctions::subRatio (type ratio1 , type 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*/
/* */
ret  ::math::rationalfunctions::multRatio (type ratio1 , type 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*/
/* */
ret  ::math::rationalfunctions::divRatio (type ratio1 , type 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 } {
 tcl = _precision 17

 f1 =     [::math::rationalfunctions::rationalFunction {1 2 3} {1 4}]
 f2 =     [::math::rationalfunctions::rationalFunction {1 2 3 0} {1 4}]
 f3 =     [::math::rationalfunctions::rationalFunction {0 0 0 0} {1}]
 f4 =     [::math::rationalfunctions::rationalFunction {5 7} {1}]
 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]"

 f1 =     [::math::rationalfunctions::rationalFunction {1 2 3} 1]
 f2 =     [::math::rationalfunctions::rationalFunction {0 1} 1]

puts "Divide:    [::math::rationalfunctions::divRatio $f1 $f2]"

 f1 =     [::math::rationalfunctions::rationalFunction {1 2 3} 1]
 f2 =     [::math::rationalfunctions::rationalFunction {1 1} {1 2}]

puts "Divide:    [::math::rationalfunctions::divRatio $f1 $f2]"

 f1 =  [::math::rationalfunctions::rationalFunction {1 2 3} 1]
 f2 =  [::math::rationalfunctions::rationalFunction {0 1} {0 0 1}]
 f3 =  [::math::rationalfunctions::divRatio $f2 $f1]
 coeffs =  [::math::rationalfunctions::coeffsNumerator $f3]
puts "Coefficients: $coeffs"
 f3 =  [::math::rationalfunctions::divRatio $f1 $f2]
 coeffs =  [::math::rationalfunctions::coeffsNumerator $f3]
puts "Coefficients: $coeffs"
 f1 =  [::math::rationalfunctions::rationalFunction {1 2 3} {1 2}]
 f2 =  [::math::rationalfunctions::rationalFunction {0} {1}]
 f3 =  [::math::rationalfunctions::divRatio $f2 $f1]
 coeffs =  [::math::rationalfunctions::coeffsNumerator $f3]
puts "Coefficients: $coeffs"
puts "Eval null function: [::math::rationalfunctions::evalRatio $f2 1]"
}