#  Created By    : Dr. Detlef Groth
 #  Created       : Mon Mar 31 09:50:42 2008
 #  Last Modified : <080403.2128>
 #
 # based on code of http://www.perlmonks.org/?node_id=482919 tlm
 # and Alan Miller's FORTRAN-implementation http://lib.stat.cmu.edu/apstat/245
 #
 ##############################################################################

 namespace eval fisher {
    set Tolerance 1
    while {[expr {1 + $Tolerance/2.0}] > 1} {
        set Tolerance [expr {$Tolerance / 2.0}]
     }
  }
 proc fisher::test {a b c d {alternative greater}} {
    if {$alternative eq "less"} {
        return [fisherExactTest $c $d $a $b]
    } else {
        return [fisherExactTest $a $b $c $d]
    }
 }
 proc fisher::fisherExactTest {a b c d {ts false}} { 
    set test [expr $a*( $a + $b + $c + $d ) - ( $a + $b )*( $a + $c )]
    if {$test < 0 && $ts} {
        return 1
    }
    if { $test < 0 } {
        if { $d < $a } {
            set p_val [_fishers_exact  $d $c $b $a $ts true]
        } else {
            set p_val [_fishers_exact  $a $b $c $d $ts true]
        }
    } else {
        if { $b < $c } {
            set p_val [_fishers_exact $b $a $d $c $ts false]
        } else {
            set p_val [_fishers_exact $c $d $a $b $ts false]
        }
    }
    return $p_val;
 }

 proc fisher::_fishers_exact {a b c d ts complement} {
    foreach {aa bb cc dd} [list $a $b $c $d] { break }
    set delta [_single_term  $aa $bb $cc $dd]
    set first $delta
    set p_val 0
    while {true} {
        set p_val [expr {$p_val + $delta}]
        if {$aa < 1} {
            break
        }
        incr bb
        incr cc
        set delta [expr {$delta * ( 1.0* ( $aa * $dd )/( $bb * $cc ) )}]
        incr aa -1
        incr dd -1
    }
    if { $ts } {
        if {$b < $c} {
            set m $b
        } else {
            set m $c
        }
        set aa [expr $a+$m]
        set bb [expr $b-$m]
        set cc [expr $c-$m]
        set dd [expr $d+$m]
        set delta  [_single_term  $aa $bb $cc $dd]
        set bound [expr {-1.0*$Tolerance}]
        while { $bound <= [expr {( $first - $delta )/$first}] && $aa > $a } {
            set p_val [expr {$p_val + $delta}];
            incr bb
            incr cc
            set delta [expr {$delta * ( ( $aa * $dd )/( $bb * $cc ) )}]
            incr aa -1
            incr dd -1
            
        }
    } elseif { $complement } {
        set p_val [expr {1 - $p_val + $first}]
    }
    return $p_val
 }
 proc fisher::_single_term {a b c d} {
    set r1 [expr {$a + $b}]
    set r2 [expr {$c + $d}]
    set c1 [expr {$a + $c}]
    set c2 [expr {$b + $d}]
    set N [expr {$r1 + $r2}]
    set s1 [expr {[_ln_fact $r1]+[_ln_fact $r2]+[_ln_fact $c1]+[_ln_fact $c2]-[_ln_fact $N]}]
    set s2 [expr {[_ln_fact $a]+[_ln_fact $b]+[_ln_fact $c]+[_ln_fact $d]}]
    return [expr {exp($s1-$s2)}]
 }
 proc fisher::_ln_fact {n} {
    if {$n <=1} { return 0 }
    return [_ln_gamm [incr n]]
 }

 proc fisher::_ln_gamm {z} {
    # Reference: "Lanczos, C. 'A precision approximation
    # of the gamma function', J. SIAM Numer. Anal., B, 1, 86-96, 1964."
    #  Translation of  Alan Miller's FORTRAN-implementation
    # See http://lib.stat.cmu.edu/apstat/245
    set x  0
    set x [expr {$x + 0.1659470187408462e-06/($z+7)}]
    set x [expr {$x + 0.9934937113930748e-05/($z+6)}]
    set x [expr {$x - 0.1385710331296526    /($z+5)}]
    set x [expr {$x + 12.50734324009056     /($z+4)}]
    set x [expr {$x - 176.6150291498386     /($z+3)}]
    set x [expr {$x + 771.3234287757674     /($z+2)}]
    set x [expr {$x - 1259.139216722289     /($z+1)}]
    set x [expr {$x + 676.5203681218835     /($z)}]
    set x [expr { $x + 0.9999999999995183}]
    return [expr {log($x)-5.58106146679532777-$z+($z-0.5)*log($z+6.5)}]
  }