Posted to tcl by kbk at Sun May 17 23:17:15 GMT 2009view raw

  1. package require math
  2. package require math::calculus
  3. for {set x 0.5} {$x <= 7.0} {set x [expr {$x + 0.5}]} {
  4. set g1 [math::calculus::romberg \
  5. [list apply {{x t} {expr {$t ** ($x-1) * exp(-$t)}}} $x] \
  6. 0 1]
  7. set g2 [math::calculus::romberg_infinity \
  8. [list apply {{x t} {expr {$t ** ($x-1) * exp(-$t)}}} $x] \
  9. 1 Inf]
  10. set gamma [expr {[lindex $g1 0] + [lindex $g2 0]}]
  11. set libgamma [expr {exp([math::ln_Gamma $x])}]
  12. puts [list $x $gamma $libgamma]
  13. }