#! perl -slw use strict; use Benchmark::Timer; use List::Util qw[ sum reduce ]; use Math::Pari qw[ factorial ]; $a=$b; sub product{ reduce{ $a *= $b } 1, @_ } sub FishersExactTest { my @data = @_; return unless @data == 4; my @C = ( sum( @data[ 0, 2 ] ), sum( @data[ 1, 3 ] ) ); my @R = ( sum( @data[ 0, 1 ] ), sum( @data[ 2, 3 ] ) ); my $N = sum @C; my $dividend = product map{ factorial $_ } grep $_, @R, @C; my $divisor = product map{ factorial $_ } grep $_, $N, @data; return $dividend / $divisor; } my $T = new Benchmark::Timer; $T->start( '' ); print FishersExactTest 989, 9400, 43300, 2400;; $T->stop( '' ); $T->report; __END__ P:\test>MP-FET.pl 8.070604647867604097576877668E-7030 1 trial of _default ( 25.852ms total), 25.852ms/trial #### #! perl -slw use strict; use Benchmark::Timer; use List::Util qw[ reduce ]; use Math::BigFloat; $a=$b; sub product{ reduce{ $a *= $b } 1, @_ } sub sum{ reduce{ $a += $b } 0, @_ } sub FishersExactTest { my @data = map{ Math::BigFloat->new( $_ ) } @_; return unless @data == 4; my @C = ( sum( @data[ 0, 2 ] ), sum( @data[ 1, 3 ] ) ); my @R = ( sum( @data[ 0, 1 ] ), sum( @data[ 2, 3 ] ) ); my $N = sum @C; my $dividend = product map{ $_->bfac } grep $_, @R, @C; my $divisor = product map{ $_->bfac } grep $_, $N, @data; return $dividend / $divisor; } my $T = new Benchmark::Timer; $T->start( '' ); print FishersExactTest 989, 9400, 43300, 2400;; $T->stop( '' ); $T->report; #### sub FishersExactTest2 { my @data = @_; return unless @data == 4; my @C = ( sum( @data[ 0, 2 ] ), sum( @data[ 1, 3 ] ) ); my @R = ( sum( @data[ 0, 1 ] ), sum( @data[ 2, 3 ] ) ); my $N = sum @C; my %dividends; $dividends{ $_ }++ for map{ factors $_ } grep $_, @R, @C; my %divisors; $divisors { $_ }++ for map{ factors $_ } grep $_, $N, @data; for my $i ( keys %divisors ) { if( exists $dividends{ $i } ) { $divisors{ $i }--, $dividends{ $i }-- while $divisors{ $i } and $dividends{ $i }; delete $divisors { $i } unless $divisors { $i }; delete $dividends{ $i } unless $dividends{ $i }; } } my $dividend = product( map{ ( $_ ) x $dividends{ $_ } } keys %dividends ); my $divisor = product( map{ ( $_ ) x $divisors { $_ } } keys %divisors ); return $dividend / $divisor; }