#! 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;
}