1# polynomials.tcl -- 2# Implement procedures to deal with polynomial functions 3# 4namespace eval ::math::polynomials { 5 variable count 0 ;# Count the number of specific commands 6 namespace eval v {} 7 8 namespace export polynomial polynCmd evalPolyn \ 9 degreePolyn coeffPolyn allCoeffsPolyn \ 10 derivPolyn primitivePolyn \ 11 addPolyn subPolyn multPolyn \ 12 divPolyn remainderPolyn 13} 14 15 16# polynomial -- 17# Return a polynomial definition 18# 19# Arguments: 20# coeffs The coefficients of the polynomial 21# Result: 22# Polynomial definition 23# 24proc ::math::polynomials::polynomial {coeffs} { 25 26 set rev_coeffs {} 27 set degree -1 28 set index 0 29 foreach coeff $coeffs { 30 if { ! [string is double -strict $coeff] } { 31 return -code error "Coefficients must be real numbers" 32 } 33 set rev_coeffs [concat $coeff $rev_coeffs] 34 if { $coeff != 0.0 } { 35 set degree $index 36 } 37 incr index 38 } 39 40 # 41 # The leading coefficient must be non-zero 42 # 43 return [list POLYNOMIAL [lrange $rev_coeffs end-$degree end]] 44} 45 46# polynCmd -- 47# Return a procedure that implements a polynomial evaluation 48# 49# Arguments: 50# coeffs The coefficients of the polynomial (or a definition) 51# Result: 52# New procedure 53# 54proc ::math::polynomials::polynCmd {coeffs} { 55 variable count 56 57 if { [lindex $coeffs 0] == "POLYNOMIAL" } { 58 set coeffs [allCoeffsPolyn $coeffs] 59 } 60 61 set degree [expr {[llength $coeffs]-1}] 62 set body "expr \{[join $coeffs +\$x*(][string repeat ) $degree]\}" 63 64 incr count 65 set name "::math::polynomials::v::POLYN$count" 66 proc $name {x} $body 67 return $name 68} 69 70# evalPolyn -- 71# Evaluate a polynomial at a given coordinate 72# 73# Arguments: 74# polyn Polynomial definition 75# x Coordinate 76# Result: 77# Value at x 78# 79proc ::math::polynomials::evalPolyn {polyn x} { 80 if { [lindex $polyn 0] != "POLYNOMIAL" } { 81 return -code error "Not a polynomial" 82 } 83 if { ! [string is double $x] } { 84 return -code error "Coordinate must be a real number" 85 } 86 87 set result 0.0 88 foreach c [lindex $polyn 1] { 89 set result [expr {$result*$x+$c}] 90 } 91 return $result 92} 93 94# degreePolyn -- 95# Return the degree of the polynomial 96# 97# Arguments: 98# polyn Polynomial definition 99# Result: 100# The degree 101# 102proc ::math::polynomials::degreePolyn {polyn} { 103 if { [lindex $polyn 0] != "POLYNOMIAL" } { 104 return -code error "Not a polynomial" 105 } 106 return [expr {[llength [lindex $polyn 1]]-1}] 107} 108 109# coeffPolyn -- 110# Return the coefficient of the index'th degree of the polynomial 111# 112# Arguments: 113# polyn Polynomial definition 114# index Degree for which to return the coefficient 115# Result: 116# The coefficient of degree "index" 117# 118proc ::math::polynomials::coeffPolyn {polyn index} { 119 if { [lindex $polyn 0] != "POLYNOMIAL" } { 120 return -code error "Not a polynomial" 121 } 122 set coeffs [lindex $polyn 1] 123 if { $index < 0 || $index > [llength $coeffs] } { 124 return -code error "Index must be between 0 and [llength $coeffs]" 125 } 126 return [lindex $coeffs end-$index] 127} 128 129# allCoeffsPolyn -- 130# Return the coefficients of the polynomial 131# 132# Arguments: 133# polyn Polynomial definition 134# Result: 135# The coefficients in ascending order 136# 137proc ::math::polynomials::allCoeffsPolyn {polyn} { 138 if { [lindex $polyn 0] != "POLYNOMIAL" } { 139 return -code error "Not a polynomial" 140 } 141 set rev_coeffs [lindex $polyn 1] 142 set coeffs {} 143 foreach c $rev_coeffs { 144 set coeffs [concat $c $coeffs] 145 } 146 return $coeffs 147} 148 149# derivPolyn -- 150# Return the derivative of the polynomial 151# 152# Arguments: 153# polyn Polynomial definition 154# Result: 155# The new polynomial 156# 157proc ::math::polynomials::derivPolyn {polyn} { 158 if { [lindex $polyn 0] != "POLYNOMIAL" } { 159 return -code error "Not a polynomial" 160 } 161 set coeffs [lindex $polyn 1] 162 set new_coeffs {} 163 set idx [degreePolyn $polyn] 164 foreach c [lrange $coeffs 0 end-1] { 165 lappend new_coeffs [expr {$idx*$c}] 166 incr idx -1 167 } 168 return [list POLYNOMIAL $new_coeffs] 169} 170 171# primitivePolyn -- 172# Return the primitive of the polynomial 173# 174# Arguments: 175# polyn Polynomial definition 176# Result: 177# The new polynomial 178# 179proc ::math::polynomials::primitivePolyn {polyn} { 180 if { [lindex $polyn 0] != "POLYNOMIAL" } { 181 return -code error "Not a polynomial" 182 } 183 set coeffs [lindex $polyn 1] 184 set new_coeffs {} 185 set idx [llength $coeffs] 186 foreach c [lrange $coeffs 0 end] { 187 lappend new_coeffs [expr {$c/double($idx)}] 188 incr idx -1 189 } 190 return [list POLYNOMIAL [concat $new_coeffs 0.0]] 191} 192 193# addPolyn -- 194# Add two polynomials and return the result 195# 196# Arguments: 197# polyn1 First polynomial or a scalar 198# polyn2 Second polynomial or a scalar 199# Result: 200# The sum of the two polynomials 201# Note: 202# Make sure that the first coefficient is not zero 203# 204proc ::math::polynomials::addPolyn {polyn1 polyn2} { 205 if { [llength $polyn1] == 1 && [string is double -strict $polyn1] } { 206 set polyn1 [polynomial $polyn1] 207 } 208 if { [llength $polyn2] == 1 && [string is double -strict $polyn2] } { 209 set polyn2 [polynomial $polyn2] 210 } 211 if { [lindex $polyn1 0] != "POLYNOMIAL" || 212 [lindex $polyn2 0] != "POLYNOMIAL" } { 213 return -code error "Both arguments must be polynomials or a real number" 214 } 215 set coeffs1 [lindex $polyn1 1] 216 set coeffs2 [lindex $polyn2 1] 217 218 set extra1 [expr {[llength $coeffs2]-[llength $coeffs1]}] 219 while { $extra1 > 0 } { 220 set coeffs1 [concat 0.0 $coeffs1] 221 incr extra1 -1 222 } 223 224 set extra2 [expr {[llength $coeffs1]-[llength $coeffs2]}] 225 while { $extra2 > 0 } { 226 set coeffs2 [concat 0.0 $coeffs2] 227 incr extra2 -1 228 } 229 230 set new_coeffs {} 231 foreach c1 $coeffs1 c2 $coeffs2 { 232 lappend new_coeffs [expr {$c1+$c2}] 233 } 234 while { [lindex $new_coeffs 0] == 0.0 } { 235 set new_coeffs [lrange $new_coeffs 1 end] 236 } 237 return [list POLYNOMIAL $new_coeffs] 238} 239 240# subPolyn -- 241# Subtract two polynomials and return the result 242# 243# Arguments: 244# polyn1 First polynomial or a scalar 245# polyn2 Second polynomial or a scalar 246# Result: 247# The difference of the two polynomials 248# Note: 249# Make sure that the first coefficient is not zero 250# 251proc ::math::polynomials::subPolyn {polyn1 polyn2} { 252 if { [llength $polyn1] == 1 && [string is double -strict $polyn1] } { 253 set polyn1 [polynomial $polyn1] 254 } 255 if { [llength $polyn2] == 1 && [string is double -strict $polyn2] } { 256 set polyn2 [polynomial $polyn2] 257 } 258 if { [lindex $polyn1 0] != "POLYNOMIAL" || 259 [lindex $polyn2 0] != "POLYNOMIAL" } { 260 return -code error "Both arguments must be polynomials or a real number" 261 } 262 set coeffs1 [lindex $polyn1 1] 263 set coeffs2 [lindex $polyn2 1] 264 265 set extra1 [expr {[llength $coeffs2]-[llength $coeffs1]}] 266 while { $extra1 > 0 } { 267 set coeffs1 [concat 0.0 $coeffs1] 268 incr extra1 -1 269 } 270 271 set extra2 [expr {[llength $coeffs1]-[llength $coeffs2]}] 272 while { $extra2 > 0 } { 273 set coeffs2 [concat 0.0 $coeffs2] 274 incr extra2 -1 275 } 276 277 set new_coeffs {} 278 foreach c1 $coeffs1 c2 $coeffs2 { 279 lappend new_coeffs [expr {$c1-$c2}] 280 } 281 while { [lindex $new_coeffs 0] == 0.0 } { 282 set new_coeffs [lrange $new_coeffs 1 end] 283 } 284 return [list POLYNOMIAL $new_coeffs] 285} 286 287# multPolyn -- 288# Multiply two polynomials and return the result 289# 290# Arguments: 291# polyn1 First polynomial or a scalar 292# polyn2 Second polynomial or a scalar 293# Result: 294# The difference of the two polynomials 295# Note: 296# Make sure that the first coefficient is not zero 297# 298proc ::math::polynomials::multPolyn {polyn1 polyn2} { 299 if { [llength $polyn1] == 1 && [string is double -strict $polyn1] } { 300 set polyn1 [polynomial $polyn1] 301 } 302 if { [llength $polyn2] == 1 && [string is double -strict $polyn2] } { 303 set polyn2 [polynomial $polyn2] 304 } 305 if { [lindex $polyn1 0] != "POLYNOMIAL" || 306 [lindex $polyn2 0] != "POLYNOMIAL" } { 307 return -code error "Both arguments must be polynomials or a real number" 308 } 309 310 set coeffs1 [lindex $polyn1 1] 311 set coeffs2 [lindex $polyn2 1] 312 313 # 314 # Take care of the null polynomial 315 # 316 if { $coeffs1 == {} || $coeffs2 == {} } { 317 return [polynomial {}] 318 } 319 320 set zeros {} 321 foreach c $coeffs1 { 322 lappend zeros 0.0 323 } 324 325 set new_coeffs [lrange $zeros 1 end] 326 foreach c $coeffs2 { 327 lappend new_coeffs 0.0 328 } 329 330 set idx 0 331 foreach c $coeffs1 { 332 set term_coeffs {} 333 foreach c2 $coeffs2 { 334 lappend term_coeffs [expr {$c*$c2}] 335 } 336 set term_coeffs [concat [lrange $zeros 0 [expr {$idx-1}]] \ 337 $term_coeffs \ 338 [lrange $zeros [expr {$idx+1}] end]] 339 340 set sum_coeffs {} 341 foreach t $term_coeffs n $new_coeffs { 342 lappend sum_coeffs [expr {$t+$n}] 343 } 344 set new_coeffs $sum_coeffs 345 incr idx 346 } 347 348 return [list POLYNOMIAL $new_coeffs] 349} 350 351# divPolyn -- 352# Divide two polynomials and return the quotient 353# 354# Arguments: 355# polyn1 First polynomial or a scalar 356# polyn2 Second polynomial or a scalar 357# Result: 358# The difference of the two polynomials 359# Note: 360# Make sure that the first coefficient is not zero 361# 362proc ::math::polynomials::divPolyn {polyn1 polyn2} { 363 if { [llength $polyn1] == 1 && [string is double -strict $polyn1] } { 364 set polyn1 [polynomial $polyn1] 365 } 366 if { [llength $polyn2] == 1 && [string is double -strict $polyn2] } { 367 set polyn2 [polynomial $polyn2] 368 } 369 if { [lindex $polyn1 0] != "POLYNOMIAL" || 370 [lindex $polyn2 0] != "POLYNOMIAL" } { 371 return -code error "Both arguments must be polynomials or a real number" 372 } 373 374 set coeffs1 [lindex $polyn1 1] 375 set coeffs2 [lindex $polyn2 1] 376 377 # 378 # Take care of the null polynomial 379 # 380 if { $coeffs1 == {} } { 381 return [polynomial {}] 382 } 383 if { $coeffs2 == {} } { 384 return -code error "Denominator can not be zero" 385 } 386 387 foreach {quotient remainder} [DivRemPolyn $polyn1 $polyn2] {break} 388 return $quotient 389} 390 391# remainderPolyn -- 392# Divide two polynomials and return the remainder 393# 394# Arguments: 395# polyn1 First polynomial or a scalar 396# polyn2 Second polynomial or a scalar 397# Result: 398# The difference of the two polynomials 399# Note: 400# Make sure that the first coefficient is not zero 401# 402proc ::math::polynomials::remainderPolyn {polyn1 polyn2} { 403 if { [llength $polyn1] == 1 && [string is double -strict $polyn1] } { 404 set polyn1 [polynomial $polyn1] 405 } 406 if { [llength $polyn2] == 1 && [string is double -strict $polyn2] } { 407 set polyn2 [polynomial $polyn2] 408 } 409 if { [lindex $polyn1 0] != "POLYNOMIAL" || 410 [lindex $polyn2 0] != "POLYNOMIAL" } { 411 return -code error "Both arguments must be polynomials or a real number" 412 } 413 414 set coeffs1 [lindex $polyn1 1] 415 set coeffs2 [lindex $polyn2 1] 416 417 # 418 # Take care of the null polynomial 419 # 420 if { $coeffs1 == {} } { 421 return [polynomial {}] 422 } 423 if { $coeffs2 == {} } { 424 return -code error "Denominator can not be zero" 425 } 426 427 foreach {quotient remainder} [DivRemPolyn $polyn1 $polyn2] {break} 428 return $remainder 429} 430 431# DivRemPolyn -- 432# Divide two polynomials and return the quotient and remainder 433# 434# Arguments: 435# polyn1 First polynomial or a scalar 436# polyn2 Second polynomial or a scalar 437# Result: 438# The difference of the two polynomials 439# Note: 440# Make sure that the first coefficient is not zero 441# 442proc ::math::polynomials::DivRemPolyn {polyn1 polyn2} { 443 444 set coeffs1 [lindex $polyn1 1] 445 set coeffs2 [lindex $polyn2 1] 446 447 set steps [expr { [degreePolyn $polyn1] - [degreePolyn $polyn2] + 1 }] 448 449 # 450 # Special case: polynomial 1 has lower degree than polynomial 2 451 # 452 if { $steps <= 0 } { 453 return [list [polynomial 0.0] $polyn1] 454 } else { 455 set extra_coeffs {} 456 for { set i 1 } { $i < $steps } { incr i } { 457 lappend extra_coeffs 0.0 458 } 459 lappend extra_coeffs 1.0 460 } 461 462 set c2 [lindex $coeffs2 0] 463 set quot_coeffs {} 464 465 for { set i 0 } { $i < $steps } { incr i } { 466 set c1 [lindex $coeffs1 0] 467 set factor [expr {$c1/$c2}] 468 469 set fpolyn [multPolyn $polyn2 \ 470 [polynomial [lrange $extra_coeffs $i end]]] 471 472 set newpol [subPolyn $polyn1 [multPolyn $fpolyn $factor]] 473 474 # 475 # Due to rounding errors, a very small, parasitical 476 # term may still exist. Remove it 477 # 478 if { [degreePolyn $newpol] == [degreePolyn $polyn1] } { 479 set new_coeffs [lrange [allCoeffsPolyn $newpol] 0 end-1] 480 set newpol [polynomial $new_coeffs] 481 } 482 set polyn1 $newpol 483 set coeffs1 [lindex $polyn1 1] 484 set quot_coeffs [concat $factor $quot_coeffs] 485 } 486 set quotient [polynomial $quot_coeffs] 487 488 return [list $quotient $polyn1] 489} 490 491# 492# Announce our presence 493# 494package provide math::polynomials 1.0.1 495 496# some tests -- 497# 498if { 0 } { 499set tcl_precision 17 500 501set f1 [::math::polynomials::polynomial {1 2 3}] 502set f2 [::math::polynomials::polynomial {1 2 3 0}] 503set f3 [::math::polynomials::polynomial {0 0 0 0}] 504set f4 [::math::polynomials::polynomial {5 7}] 505set cmdf1 [::math::polynomials::polynCmd {1 2 3}] 506 507foreach x {0 1 2 3 4 5} { 508 puts "[::math::polynomials::evalPolyn $f1 $x] -- \ 509[expr {1.0+2.0*$x+3.0*$x*$x}] -- \ 510[$cmdf1 $x] -- [::math::polynomials::evalPolyn $f3 $x]" 511} 512 513puts "Degree: [::math::polynomials::degreePolyn $f1] (expected: 2)" 514puts "Degree: [::math::polynomials::degreePolyn $f2] (expected: 2)" 515foreach d {0 1 2} { 516 puts "Coefficient $d = [::math::polynomials::coeffPolyn $f2 $d]" 517} 518puts "All coefficients = [::math::polynomials::allCoeffsPolyn $f2]" 519 520puts "Derivative = [::math::polynomials::derivPolyn $f1]" 521puts "Primitive = [::math::polynomials::primitivePolyn $f1]" 522 523puts "Add: [::math::polynomials::addPolyn $f1 $f4]" 524puts "Add: [::math::polynomials::addPolyn $f4 $f1]" 525puts "Subtract: [::math::polynomials::subPolyn $f1 $f4]" 526puts "Multiply: [::math::polynomials::multPolyn $f1 $f4]" 527 528set f1 [::math::polynomials::polynomial {1 2 3}] 529set f2 [::math::polynomials::polynomial {0 1}] 530 531puts "Divide: [::math::polynomials::divPolyn $f1 $f2]" 532puts "Remainder: [::math::polynomials::remainderPolyn $f1 $f2]" 533 534set f1 [::math::polynomials::polynomial {1 2 3}] 535set f2 [::math::polynomials::polynomial {1 1}] 536 537puts "Divide: [::math::polynomials::divPolyn $f1 $f2]" 538puts "Remainder: [::math::polynomials::remainderPolyn $f1 $f2]" 539 540set f1 [::math::polynomials::polynomial {1 2 3}] 541set f2 [::math::polynomials::polynomial {0 1}] 542set f3 [::math::polynomials::divPolyn $f2 $f1] 543set coeffs [::math::polynomials::allCoeffsPolyn $f3] 544puts "Coefficients: $coeffs" 545set f3 [::math::polynomials::divPolyn $f1 $f2] 546set coeffs [::math::polynomials::allCoeffsPolyn $f3] 547puts "Coefficients: $coeffs" 548set f1 [::math::polynomials::polynomial {1 2 3}] 549set f2 [::math::polynomials::polynomial {0}] 550set f3 [::math::polynomials::divPolyn $f2 $f1] 551set coeffs [::math::polynomials::allCoeffsPolyn $f3] 552puts "Coefficients: $coeffs" 553} 554