#!/usr/bin/perl

use Math::Complex;
use Math::Polynomial::Solve qw(:classical coefficients);
coefficients order => 'descending';

open IN,"0.asc";

$four=0;
$three=0;
$repeat=0;
$maxerror = 0;
$maxerror3 = 0;
$gmaxerror = 0;
$gmaxerror3 = 0;

while($l = <IN>)
{
  chomp($l);
# print "Line: $l\n";
  if($four)
  {
#   print "four: $l\n";
    if($l =~ /^  (.{13})  (.{13})  (.{13})$/)
    {
    $real=$1;
    $imag=$2;
    $err=$3;
    $real =~ s/'/e/g;
    $real =~ s/ //g;
    $real += 0;
    $reals[$repeat][$four] = $real;
    $imag =~ s/'/e/g;
    $imag =~ s/ //g;
    $imag += 0;
    $imags[$repeat][$four] = $imag;
    $err =~ s/'/e/g;
    $err =~ s/ //g;
    $err += 0;
    $errs[$repeat][$four] = $err;
    $gmaxerror = $err if($err>$gmaxerror);
#   print "four I: $real $imag $err\n";
    }
    $four++;
    if($four==5)
    {
      $four = 0;
      for($i=0;$i<=4;$i++)
      {
	printf "COEF[%d]: %20.12e\n",$i,$coef[$repeat][$i];
      }
      print "\nSolution type: ".$soltype[$repeat]."\n";
      for($i=1;$i<=4;$i++)
      {
	printf "Real: %14.10f imag: %14.10f err: %14.10e\n",$reals[$repeat][$i],$imags[$repeat][$i],$errs[$repeat][$i];
      }
      print "\n\n";
      @x4 = quartic_roots($coef[$repeat][4], $coef[$repeat][3], $coef[$repeat][2], $coef[$repeat][1], $coef[$repeat][0]);
      $nroot = scalar @x4;
      $soltype2=0;
      for($i=0;$i<=3;$i++)
      {
	$x = $x4[$i];
	$err = 0;
	for($j=4;$j>=0;$j--)
	{
	  $err = $err*$x+$coef[$repeat][$j];
	}
	$err = abs($err);
	$maxerror = $err if($err > $maxerror);
	if(Im($x4[$i])!=0)
	{
	  $soltype2 += 1;
	}
	else
	{
	  $soltype2 += 10;
	}
	print $x4[$i]."\t".$err."\n";
      }
      if($soltype[$repeat] == $soltype2)
      {
	print "soltype match\n";
      }
      else
      {
	print "soltype MISMATCH\n";
      }
    }
    next;
  }
  if($three)
  {
    print "three $l\n";
    if($l =~ /^ \d  (.{13})  (.{13})  (.{13})$/)
    {
    $real=$1;
    $imag=$2;
    $err=$3;
    $real =~ s/'/e/g;
    $real =~ s/ //g;
    $real += 0;
    $reals3[$repeat][$three] = $real;
    $imag =~ s/'/e/g;
    $imag =~ s/ //g;
    $imag += 0;
    $imags3[$repeat][$three] = $imag;
    $err =~ s/'/e/g;
    $err =~ s/ //g;
    $err += 0;
    $errs3[$repeat][$three] = $err;
    $gmaxerror3 = $err if($err>$gmaxerror3);
    print "three I: $real $imag $err\n";
    }
    $three++;
    if($three==4)
    {
      $three = 0;
      for($i=0;$i<=3;$i++)
      {
	printf "COEF3[%d]: %20.12e\n",$i,$coef3[$repeat][$i];
      }
      for($i=1;$i<=3;$i++)
      {
	printf "Real: %14.10f imag: %14.10f\n",$reals3[$repeat][$i],$imags3[$repeat][$i];
      }
      print "\n\n";
      @x3 = cubic_roots($coef3[$repeat][3], $coef3[$repeat][2], $coef3[$repeat][1], $coef3[$repeat][0]);
      $nroot = scalar @x3;
      for($i=0;$i<=2;$i++)
      {
	$x = $x3[$i];
	$err = 0;
	for($j=3;$j>=0;$j--)
	{
	  $err = $err*$x+$coef3[$repeat][$j];
	}
	$err = abs($err);
	$maxerror3 = $err if($err > $maxerror3);
	print $x3[$i]."\t$err\n";
      }
    }
    next;
  }
  if($l =~ /^repeat (\d+)$/)
  {
    $repeat = $1+0;
    print "============== $repeat ================\n";
  }
  $c = substr($l,2,1);
  if($l =~ /^cubic (\d)  (.{14})$/)
  {
    $icoef3 = $1;
    $type2 = $2;
    $type2 =~ s/'/e/g;
    $type2 =~ s/ //g;
    $type2 += 0;
    $coef3[$repeat][$icoef3] = $type2;
    print "coef3: $icoef3 $type2\n";
  }
  if($l =~ /^ ([0-4])  (.*)$/)
  {
    $icoef = $1;
    $type2 = $2;
    $type2 =~ s/'/e/g;
    $type2 =~ s/ //g;
    $type2 += 0;
    $coef[$repeat][$icoef] = $type2;
#   print "kilroy 2: $type2\n";
  }
  if($l =~ /^branch (\d)$/)
  {
    $branchtype = $1;
    print "branchtype:  $branchtype\n";
    $three=1;
  }
  if($l =~ /^(.*) roots$/)
  {
    $nsol = $1;
    if($nsol eq "four complex")
    {
      $soltype[$repeat]=4;
    }
    elsif($nsol eq "two real and two complex")
    {
      $soltype[$repeat]=22;
    }
    elsif($nsol eq "four real")
    {
      $soltype[$repeat]=40;
    }
    else
    {
      $soltype[$repeat]=-1;
      print "Unknown: $nsol\n";
    }

#   print "kilroy 3: $nsol $soltype[$repeat]\n";
    $four = 1;
  }
}

print "maxerror: $maxerror\n";
print "maxerror3: $maxerror3\n";
print "gmaxerror: $gmaxerror\n";
print "gmaxerror3: $gmaxerror3\n";

$k=0;
for($i=1; $i<=$repeat; $i++)
{
  for($j=1; $j<=4; $j++)
  {
    $totalerr[$k] = abs($errs[$i][$j]);
    $totalno[$k] = $i;
    $totalroot[$k] = $j;
    $k++;
  }
}

@sorted_indices = sort { $totalerr[$b] <=> $totalerr[$a] } 0..$#totalerr;

for($i=0; $i<$#totalerr; $i++)
{
  $j = $sorted_indices[$i];
  print $totalerr[$j]."\t".$totalno[$j]."\t".$totalroot[$j]."\n";
}