1# rational_funcs.tcl -- 2# Implement procedures to deal with rational functions 3# 4 5package require math::polynomials 6 7namespace eval ::math::rationalfunctions { 8 variable count 0 ;# Count the number of specific commands 9 namespace eval v {} 10 11 namespace export rationalFunction ratioCmd evalRatio \ 12 coeffsNumerator coeffsDenominator \ 13 derivRatio \ 14 addRatio subRatio multRatio \ 15 divRatio 16 17 namespace import ::math::polynomials::* 18} 19 20 21# rationalFunction -- 22# Return a rational function definition 23# 24# Arguments: 25# num The coefficients of the numerator 26# den The coefficients of the denominator 27# Result: 28# Rational function definition 29# 30proc ::math::rationalfunctions::rationalFunction {num den} { 31 32 foreach coeffs [list $num $den] { 33 foreach coeff $coeffs { 34 if { ! [string is double -strict $coeff] } { 35 return -code error "Coefficients must be real numbers" 36 } 37 } 38 } 39 40 # 41 # The leading coefficient must be non-zero 42 # 43 return [list RATIONAL_FUNCTION [polynomial $num] [polynomial $den]] 44} 45 46# ratioCmd -- 47# Return a procedure that implements a rational function evaluation 48# 49# Arguments: 50# num The coefficients of the numerator 51# den The coefficients of the denominator 52# Result: 53# New procedure 54# 55proc ::math::rationalfunctions::ratioCmd {num {den {}}} { 56 variable count 57 58 if { [llength $den] == 0 } { 59 if { [lindex $num 0] == "RATIONAL_FUNCTION" } { 60 set den [lindex $num 2] 61 set num [lindex $num 1] 62 } 63 } 64 65 set degree1 [expr {[llength $num]-1}] 66 set degree2 [expr {[llength $num]-1}] 67 set body "expr \{([join $num +\$x*(][string repeat ) $degree1])/\ 68(double([join $den +\$x*(][string repeat ) $degree2])\}" 69 70 incr count 71 set name "::math::rationalfunctions::v::RATIO$count" 72 proc $name {x} $body 73 return $name 74} 75 76# evalRatio -- 77# Evaluate a rational function at a given coordinate 78# 79# Arguments: 80# ratio Rational function definition 81# x Coordinate 82# Result: 83# Value at x 84# 85proc ::math::rationalfunctions::evalRatio {ratio x} { 86 if { [lindex $ratio 0] != "RATIONAL_FUNCTION" } { 87 return -code error "Not a rational function" 88 } 89 if { ! [string is double $x] } { 90 return -code error "Coordinate must be a real number" 91 } 92 93 set num 0.0 94 foreach c [lindex [lindex $ratio 1] 1] { 95 set num [expr {$num*$x+$c}] 96 } 97 98 set den 0.0 99 foreach c [lindex [lindex $ratio 2] 1] { 100 set den [expr {$den*$x+$c}] 101 } 102 return [expr {$num/double($den)}] 103} 104 105# coeffsNumerator -- 106# Return the coefficients of the numerator 107# 108# Arguments: 109# ratio Rational function definition 110# Result: 111# The coefficients in ascending order 112# 113proc ::math::rationalfunctions::coeffsNumerator {ratio} { 114 if { [lindex $ratio 0] != "RATIONAL_FUNCTION" } { 115 return -code error "Not a rational function" 116 } 117 set polyn [lindex $ratio 1] 118 return [allCoeffsPolyn $polyn] 119} 120 121# coeffsDenominator -- 122# Return the coefficients of the denominator 123# 124# Arguments: 125# ratio Rational function definition 126# Result: 127# The coefficients in ascending order 128# 129proc ::math::rationalfunctions::coeffsDenominator {ratio} { 130 if { [lindex $ratio 0] != "RATIONAL_FUNCTION" } { 131 return -code error "Not a rational function" 132 } 133 set polyn [lindex $ratio 2] 134 return [allCoeffsPolyn $polyn] 135} 136 137# derivRatio -- 138# Return the derivative of the rational function 139# 140# Arguments: 141# polyn Polynomial definition 142# Result: 143# The new polynomial 144# 145proc ::math::rationalfunctions::derivRatio {ratio} { 146 if { [lindex $ratio 0] != "RATIONAL_FUNCTION" } { 147 return -code error "Not a rational function" 148 } 149 set num_polyn [lindex $ratio 1] 150 set den_polyn [lindex $ratio 2] 151 set num_deriv [derivPolyn $num_polyn] 152 set den_deriv [derivPolyn $den_polyn] 153 set num [subPolyn [multPolyn $num_deriv $den_polyn] \ 154 [multPolyn $den_deriv $num_polyn] ] 155 set den [multPolyn $den_polyn $den_polyn] 156 157 return [list RATIONAL_FUNCTION $num $den] 158} 159 160# addRatio -- 161# Add two rational functions and return the result 162# 163# Arguments: 164# ratio1 First rational function or a scalar 165# ratio2 Second rational function or a scalar 166# Result: 167# The sum of the two functions 168# Note: 169# TODO: Check for the same denominator 170# 171proc ::math::rationalfunctions::addRatio {ratio1 ratio2} { 172 if { [llength $ratio1] == 1 && [string is double -strict $ratio1] } { 173 set polyn1 [rationalFunction $ratio1 1.0] 174 } 175 if { [llength $ratio2] == 1 && [string is double -strict $ratio2] } { 176 set ratio2 [rationalFunction $ratio1 1.0] 177 } 178 if { [lindex $ratio1 0] != "RATIONAL_FUNCTION" || 179 [lindex $ratio2 0] != "RATIONAL_FUNCTION" } { 180 return -code error "Both arguments must be rational functions or a real number" 181 } 182 183 set num1 [lindex $ratio1 1] 184 set den1 [lindex $ratio1 2] 185 set num2 [lindex $ratio2 1] 186 set den2 [lindex $ratio2 2] 187 188 set newnum [addPolyn [multPolyn $num1 $den2] \ 189 [multPolyn $num2 $den1] ] 190 191 set newden [multPolyn $den1 $den2] 192 193 return [list RATIONAL_FUNCTION $newnum $newden] 194} 195 196# subRatio -- 197# Subtract two rational functions and return the result 198# 199# Arguments: 200# ratio1 First rational function or a scalar 201# ratio2 Second rational function or a scalar 202# Result: 203# The difference of the two functions 204# Note: 205# TODO: Check for the same denominator 206# 207proc ::math::rationalfunctions::subRatio {ratio1 ratio2} { 208 if { [llength $ratio1] == 1 && [string is double -strict $ratio1] } { 209 set polyn1 [rationalFunction $ratio1 1.0] 210 } 211 if { [llength $ratio2] == 1 && [string is double -strict $ratio2] } { 212 set ratio2 [rationalFunction $ratio1 1.0] 213 } 214 if { [lindex $ratio1 0] != "RATIONAL_FUNCTION" || 215 [lindex $ratio2 0] != "RATIONAL_FUNCTION" } { 216 return -code error "Both arguments must be rational functions or a real number" 217 } 218 219 set num1 [lindex $ratio1 1] 220 set den1 [lindex $ratio1 2] 221 set num2 [lindex $ratio2 1] 222 set den2 [lindex $ratio2 2] 223 224 set newnum [subPolyn [multPolyn $num1 $den2] \ 225 [multPolyn $num2 $den1] ] 226 227 set newden [multPolyn $den1 $den2] 228 229 return [list RATIONAL_FUNCTION $newnum $newden] 230} 231 232# multRatio -- 233# Multiply two rational functions and return the result 234# 235# Arguments: 236# ratio1 First rational function or a scalar 237# ratio2 Second rational function or a scalar 238# Result: 239# The product of the two functions 240# Note: 241# TODO: Check for the same denominator 242# 243proc ::math::rationalfunctions::multRatio {ratio1 ratio2} { 244 if { [llength $ratio1] == 1 && [string is double -strict $ratio1] } { 245 set polyn1 [rationalFunction $ratio1 1.0] 246 } 247 if { [llength $ratio2] == 1 && [string is double -strict $ratio2] } { 248 set ratio2 [rationalFunction $ratio1 1.0] 249 } 250 if { [lindex $ratio1 0] != "RATIONAL_FUNCTION" || 251 [lindex $ratio2 0] != "RATIONAL_FUNCTION" } { 252 return -code error "Both arguments must be rational functions or a real number" 253 } 254 255 set num1 [lindex $ratio1 1] 256 set den1 [lindex $ratio1 2] 257 set num2 [lindex $ratio2 1] 258 set den2 [lindex $ratio2 2] 259 260 set newnum [multPolyn $num1 $num2] 261 set newden [multPolyn $den1 $den2] 262 263 return [list RATIONAL_FUNCTION $newnum $newden] 264} 265 266# divRatio -- 267# Divide two rational functions and return the result 268# 269# Arguments: 270# ratio1 First rational function or a scalar 271# ratio2 Second rational function or a scalar 272# Result: 273# The quotient of the two functions 274# Note: 275# TODO: Check for the same denominator 276# 277proc ::math::rationalfunctions::divRatio {ratio1 ratio2} { 278 if { [llength $ratio1] == 1 && [string is double -strict $ratio1] } { 279 set polyn1 [rationalFunction $ratio1 1.0] 280 } 281 if { [llength $ratio2] == 1 && [string is double -strict $ratio2] } { 282 set ratio2 [rationalFunction $ratio1 1.0] 283 } 284 if { [lindex $ratio1 0] != "RATIONAL_FUNCTION" || 285 [lindex $ratio2 0] != "RATIONAL_FUNCTION" } { 286 return -code error "Both arguments must be rational functions or a real number" 287 } 288 289 set num1 [lindex $ratio1 1] 290 set den1 [lindex $ratio1 2] 291 set num2 [lindex $ratio2 1] 292 set den2 [lindex $ratio2 2] 293 294 set newnum [multPolyn $num1 $den2] 295 set newden [multPolyn $num2 $den1] 296 297 return [list RATIONAL_FUNCTION $newnum $newden] 298} 299 300# 301# Announce our presence 302# 303package provide math::rationalfunctions 1.0.1 304 305# some tests -- 306# 307if { 0 } { 308set prec $::tcl_precision 309if {![package vsatisfies [package provide Tcl] 8.5]} { 310 set ::tcl_precision 17 311} else { 312 set ::tcl_precision 0 313} 314 315 316set f1 [::math::rationalfunctions::rationalFunction {1 2 3} {1 4}] 317set f2 [::math::rationalfunctions::rationalFunction {1 2 3 0} {1 4}] 318set f3 [::math::rationalfunctions::rationalFunction {0 0 0 0} {1}] 319set f4 [::math::rationalfunctions::rationalFunction {5 7} {1}] 320set cmdf1 [::math::rationalfunctions::ratioCmd {1 2 3} {1 4}] 321 322foreach x {0 1 2 3 4 5} { 323 puts "[::math::rationalfunctions::evalRatio $f1 $x] -- \ 324[expr {(1.0+2.0*$x+3.0*$x*$x)/double(1.0+4.0*$x)}] -- \ 325[$cmdf1 $x] -- [::math::rationalfunctions::evalRatio $f3 $x]" 326} 327 328puts "All coefficients = [::math::rationalfunctions::coeffsNumerator $f2]" 329puts " [::math::rationalfunctions::coeffsDenominator $f2]" 330 331puts "Derivative = [::math::rationalfunctions::derivRatio $f1]" 332 333puts "Add: [::math::rationalfunctions::addRatio $f1 $f4]" 334puts "Add: [::math::rationalfunctions::addRatio $f4 $f1]" 335puts "Subtract: [::math::rationalfunctions::subRatio $f1 $f4]" 336puts "Multiply: [::math::rationalfunctions::multRatio $f1 $f4]" 337 338set f1 [::math::rationalfunctions::rationalFunction {1 2 3} 1] 339set f2 [::math::rationalfunctions::rationalFunction {0 1} 1] 340 341puts "Divide: [::math::rationalfunctions::divRatio $f1 $f2]" 342 343set f1 [::math::rationalfunctions::rationalFunction {1 2 3} 1] 344set f2 [::math::rationalfunctions::rationalFunction {1 1} {1 2}] 345 346puts "Divide: [::math::rationalfunctions::divRatio $f1 $f2]" 347 348set f1 [::math::rationalfunctions::rationalFunction {1 2 3} 1] 349set f2 [::math::rationalfunctions::rationalFunction {0 1} {0 0 1}] 350set f3 [::math::rationalfunctions::divRatio $f2 $f1] 351set coeffs [::math::rationalfunctions::coeffsNumerator $f3] 352puts "Coefficients: $coeffs" 353set f3 [::math::rationalfunctions::divRatio $f1 $f2] 354set coeffs [::math::rationalfunctions::coeffsNumerator $f3] 355puts "Coefficients: $coeffs" 356set f1 [::math::rationalfunctions::rationalFunction {1 2 3} {1 2}] 357set f2 [::math::rationalfunctions::rationalFunction {0} {1}] 358set f3 [::math::rationalfunctions::divRatio $f2 $f1] 359set coeffs [::math::rationalfunctions::coeffsNumerator $f3] 360puts "Coefficients: $coeffs" 361puts "Eval null function: [::math::rationalfunctions::evalRatio $f2 1]" 362 363set ::tcl_precision $prec 364} 365