I never doubted the idea, nor your proof, I simply had real trouble coding it. There's a load of niggly edge cases that I could not seem to get right--but now I have.
As a result, the original testcase that took M::BF 4 1/2 hours and that M::Pari couldn't handle, I now have down to 192 milliseconds in pure Perl! What's more, it happily handles a dataset of (1e6 2e6 2e6 1e6) without blowing the memory in just over 1 minute (though I cannot check the accuracy as I have nothing else that will touch those numbers).
Of course, I've given up a fair amount of precision along the way, so now I am wondering if I can enhance my cheap BigFloat code to recover some of it?
#! perl -slw
use strict;
use Benchmark::Timer; my $T = new Benchmark::Timer;
use List::Util qw[ sum reduce max ]; our( $a, $b );
sub factors{ 2 .. $_[ 0 ] }
sub normalise {
my( $s, $n ) = @{+shift };
while( 1 ) {
if( $n > 1.0 ) { $n /= 10; $s++; redo; }
elsif( $n < 0.1 ) { $n *= 10; $s--; redo; }
else { last }
}
return [ $s, $n ];
}
sub sProduct{
@{
reduce{
$a->[ 0 ] += $b->[ 0 ];
$a->[ 1 ] *= $b->[ 1 ];
normalise( $a );
} @_
}
}
sub sProdRange {
map{
reduce{
$a->[ 1 ] *= $b;
normalise( $a );
} [ 0, 1 ], $_->[ 0 ] .. $_->[ 1 ]
} @_;
}
sub reduceRanges{
our( @v, @s );
local( *v, *s ) = @_;
@v = sort{ $b <=> $a } @v;
@s = sort{ $b <=> $a } @s;
for my $i ( 0 .. max( $#v, $#s ) ) {
$s[$i]=1 unless defined $s[$i];
$v[$i]=1 unless defined $v[$i];
my $cmp = $v[$i] <=> $s[$i];
unless( $cmp ) {
$s[$i] = [ 1, 1 ];
$v[$i] = [ 1, 1 ];
}
elsif( $cmp < 0 ) {
$s[$i] = [ $s[$i] - ( $s[$i] - 1 - $v[$i] ), $s[$i] ];
$v[$i] = [ 1, 1 ];
}
else {
$v[$i] = [ $v[$i] - ( $v[$i] - 1 - $s[$i] ), $v[$i] ];
$s[$i] = [ 1, 1 ];
}
}
return \@v, \@s;
}
sub FET6 {
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( $Vref, $Sref ) = reduceRanges [grep $_, @R, @C], [grep $_, $N,
+ @data];
my( $dScale, $d ) = sProduct sProdRange @$Vref;
my( $sScale, $s ) = sProduct sProdRange @$Sref;
return ( $d / $s, $dScale - $sScale );
}
die "Bad args @ARGV" unless @ARGV == 4;
$T->start("[@ARGV]");
printf "%.17fe%d\n", FET6 @ARGV;
$T->stop("[@ARGV]");
$T->report;
__END__
[ 2:39:39.81] P:\test>FET6 5 0 1 4
0.23809523809523819e-1
1 trial of [5 0 1 4] ( 728us total), 728us/trial
[ 2:39:32.26] P:\test>fet6 989 9400 43300 2400
0.80706046478677251e-7029
1 trial of [989 9400 43300 2400] (192.499ms total), 192.499ms/trial
[ 2:37:08.39] P:\test>FET6 1e6 2e6 2e6 1e6
2.73353759185676100e-147576
1 trial of [1e6 2e6 2e6 1e6] ( 67.960s total), 67.960s/trial
[ 2:39:18.10] P:\test>FET6 1e5 2e5 2e5 1e5
1.32499459649680040e-14760
1 trial of [1e5 2e5 2e5 1e5] ( 6.382s total), 6.382s/trial
[ 2:35:53.28] P:\test>FET6 1e6 2e6 2e6 1e6
[1e6 2e6 2e6 1e6]
2.73353759185676100e-147576
1 trial of _default ( 67.887s total), 67.887s/trial
