# TAOCP Vol 1, Section 2.2.6, Exercise 15. MMIX version. 
# Perl wrapping for testing.

use  Math::MatrixReal;		# www.cpan.org

$path_to_mmix = "../../support/mmix_with_printfloat/mmix_with_printfloat";


# Sequence of random matrices.
$start = $ARGV[0] || 1;
$stop = $ARGV[1] || 5;
$dry_until = $ARGV[2] || $start;
$seed = $ARGV[3] || 0;
$step = $ARGV[4] || 1;

srand($seed);
for (my $i = $start; $i <= $stop; $i += $step) {
    my $ma = &random_matrix($i);
    if ($i >= $dry_until) {
	&run($ma);
    }
}

if (0) {
# Original matrix from the book. 
my $m1 = Math::MatrixReal->new_from_rows( [ [1,2,3], [0,1,2], [0,0,1] ] );
printf "original "; &run($m1);

# Singular matrix.
my $m2 = Math::MatrixReal->new_from_rows( [ [1,2,3], [1,2,3], [0,0,1] ] );
printf "singular matrix "; &run($m2);
}


sub random_matrix {
    my $n = shift;
    my $ma = new Math::MatrixReal($n,$n);
    for (;;) {
	for (1..$n) { 
	    &add_random_element($ma,$n);
	}
	last if abs($ma->det()) > 1e-5;
    }
    $ma;
}

sub add_random_element {
    my $ma = shift;
    my $n = shift;
    my $i = int(1 + rand($n));
    my $j = int(1 + rand($n));
    my $val = (rand() - 0.5);
    $ma->assign($i,$j,$val);
#    print "adding $i $j $val\n";
}

sub dump_matrix {
    my $ma= shift;
    my $F = shift;
    my ($m,$n) = $ma->dim();
    print $F qq{\n:inputdata	OCTA	$n\n};
    for my $i (1..$m) {
	for my $j  (1..$n) {
	    my $e = $ma->element($i,$j);
	    if ($e) {
		printf $F qq{	OCTA	%d,%d\n1H	BYTE	"%lf",0; LOC 1B+32\n},$i,$j,$e;
	    }
	}
    }
    print $F qq{	OCTA	-1\n};
}

sub create_source {
    my $ma= shift;
    my $f,$o;
    open($f,"pivot_run.mms") || die;
    open($o,">test.mms") || die;
    while(<$f>) {
	last if (/^\:inputdata/);
	print $o $_;
    }
    &dump_matrix($ma,$o);
    close($o);
    close($f);
}


sub run {
    my $ma = shift;
    my ($m,$n) = $ma->dim();
    ($m == $n) || die;

    &create_source($ma);
    system("mmixal test.mms");
    my $r,@mo;
    $mo[0] = Math::MatrixReal->new($n,$n);
    $mo[1] = Math::MatrixReal->new($n,$n);
    my $io = 0;
    open ($r,"$path_to_mmix test|") || die;
    while (<$r>) {
	if (/Input matrix/) {
	    $io= 0;
	} elsif (/Output matrix/) {
	    $io= 1;
	} elsif (/^(\d+),(\d+),([0-9.eE\-]+)/) {
	    $mo[$io]->assign($1,$2,$3);
	}
    }
    close($r);
    my $unary = Math::MatrixReal->new($m,$n)->each( sub { $_[1] == $_[2];  } );
    $err = abs($mo[0]*$mo[1]-$unary);

    printf "n=%4d err=%9.3g\n",$n,$err;

}