#!/usr/bin/perl

# double-torus-relations: given a set of t/i numbers
# ($t/$i, $u/$j, $v/$k), first map the corresponding
# knot into a path of handle/row/column numbers
# ($h, $r, $c), and then to a set of relations.

# Whether or not to output verbose intermediate step
# information for debugging
$DEBUG = 0;

use CGI;

open STDERR, '>&STDOUT';

my $q=CGI->new;

$t1 = $q->param('t1');
$i1 = $q->param('i1');
$t2 = $q->param('t2');
$i2 = $q->param('i2');
$t3 = $q->param('t3');
$i3 = $q->param('i3');

my @ti = ($t1,$i1,$t2,$i2,$t3,$i3);
print $q->header();
print '<pre>';

# el letters, for now.
my @ELS = ('a' .. 'd', 'f' .. 'z');

my @t = @ti[0, 2, 4];
my @i = @ti[1, 3, 5];

unless ( $i1 =~ /^[+-]?[0-9]+$/ ) {
  error( "i1 is not an integer.\n" );
}

unless ( $t1 =~ /^[+-]?[0-9]+$/ ) {
  error( "t1 is not an integer.\n" );
}

unless ( $i2 =~ /^[+-]?[0-9]+$/ ) {
  error( "i2 is not an integer.\n" );
}

unless ( $t2 =~ /^[+-]?[0-9]+$/ ) {
  error( "t2 is not an integer.\n" );
}

unless ( $i3 =~ /^[+-]?[0-9]+$/ ) {
  error( "i3 is not an integer.\n" );
}

unless ( $t3 =~ /^[+-]?[0-9]+$/ ) {
  error( "t3 is not an integer.\n" );
}

# immediately check if we have a 0 denom with a non-zero num.
if ($i[0] == 0 && $t[0] != 0 or
    $i[1] == 0 && $t[1] != 0 or
    $i[2] == 0 && $t[2] != 0) {
  error("a/0, with a != 0: must be not-knot.\n");
}

# verify that @i are non-negative and sum of @i is even and middle $i is maximum
verify_i(@i);

# total # of crossings.
my $crossings = abs($t[0]) * ($i[0] - 1)
              + abs($t[1]) * ($i[1] - 1)
              + abs($t[2]) * ($i[2] - 1);

# make the handle-map.
my %h_map = new_h_map( $t[0], $t[1], $t[2], $i[0], $i[1], $i[2] );
if ( $DEBUG ) 
{
    print( "handle map:\n" );
    print( h_map_str(%h_map), "\n" );
    print_map(%h_map);
    print( "\n" );
}

# make the path.
if( $DEBUG ) { print "path:\n" };

my($els, $path, $tops, %crossings) = path(\@t, \@i, \%h_map);

if( $DEBUG )
{
    print "els(@$els)\n";
    print "\n";
}

# error-checking.
{
  # XXX check this.
  if (@$els == 1 or $els->[-1] != 0) {
    my $message = "DID NOT CLOSE PATH: ELS(@$els)\n";
    error($message, $els, $path, %crossings);
  }
  pop @$els;
  if (@$els == 1) {
    # one element; must be unknot/link/unlink.
    my $message = "ONLY ONE ELEMENT\n";
    error($message, $els, $path, %crossings);
  }
  # this checks that we visited all top strands.
  if ($tops != $i[0] + $i[1] + $i[2]) {
    my $i_sum = $i[0] + $i[1] + $i[2];
    my $message = "DID NOT VISIT ALL TOP STRANDS ($tops < $i_sum)\n";
    error($message, $els, $path, %crossings);
  }
}

if( $DEBUG )
{
	print "crossings:\n";
	print_crossings(%crossings);
	print "\n";
}

# get the relators.
if( $DEBUG ) { print "relators:\n"; }
my @relators = relators(%crossings);
if( $DEBUG ) { print "\n"; }

if ( $DEBUG )
{
	print "relators again:\n";
    for my $ret (@relators)
    {
        print "$ret\n";
    }
print "\n";
}

# make the Jacobian matrix.
if ( $DEBUG ) { print "matrix:\n"; }
my @matrix = jacobian_matrix(@relators);
if( $DEBUG )
{
    print "\n";
    print "\n";
}

# finally, the real Jacobian.
print_matrix(@matrix);
print "determinant:\n";
print "\n";
my $det = det(@matrix);
$det = poly_prune($det);
$det = poly_center($det);
my $det_str = poly_str(%$det);
print "DET: $det_str\n";
my $rs = rolfsen_str($det);
print "RS: $rs\n";

exit;

sub new_h_map
{
	my ( $a, $b, $c, $x, $y, $z ) = @_;

    my $l = ( ( -( $x ) + $y - $z ) / 2 );
    
    if( $l < 0 ) { $l = 0; }
			
	my $xy = ( ( $x + $y - $z ) / 2 ) - $l;
	my $xz = ( ( $x - $y + $z ) / 2 ) + $l;
	my $yz = ( ( -( $x ) + $y + $z ) / 2) - $l;
	
	my $n = 0;
		
	my %h_map;
	
	if( $l == 0 )
	{	
		do
		{
			$h_map{ cs( [ 0, 0, $n ] ) } = cs( [ 2, 0, $z - $n - 1 ] );
			$h_map{ cs( [ 2, 0, $z - $n - 1 ] ) } = cs( [ 0, 0, $n ] );
			$h_map{ cs( [ 0, abs( $a ) + 1, $n ] ) } = cs( [ 2, abs( $c ) + 1, $z - $n - 1 ] );
			$h_map{ cs( [ 2, abs( $c ) + 1, $z - $n - 1 ] ) } = cs( [ 0, abs( $a ) + 1, $n ] );
			$n++;
		} while $n < $xz;
		
		$n = 0;
		
		do
		{
			$h_map{ cs( [ 0, 0, $xz + $n ] ) } = cs( [ 1, 0, $xy - $n - 1 ] );
			$h_map{ cs( [ 1, 0, $xy - $n - 1 ] ) } = cs( [ 0, 0, $xz + $n ] );
			$h_map{ cs( [ 0, abs( $a ) + 1, $xz + $n ] ) } = cs( [ 1, abs( $b ) + 1, $xy - $n - 1 ] );
			$h_map{ cs( [ 1, abs( $b ) + 1, $xy - $n - 1 ] ) } = cs( [ 0, abs( $a ) + 1, $xz + $n ] );
			$n++;
		} while $n < $xy;
		
		$n = 0;
		
		do
		{
			$h_map{ cs( [ 1, 0, $xy + $n ] ) } = cs( [ 2, 0, $yz - $n - 1 ] );
			$h_map{ cs( [ 2, 0, $yz - $n - 1 ] ) } = cs( [ 1, 0, $xy + $n ] );
			$h_map{ cs( [ 1, abs( $b ) + 1, $xy + $n ] ) } = cs( [ 2, abs( $c ) + 1, $yz - $n - 1 ] );
			$h_map{ cs( [ 2, abs( $c ) + 1, $yz - $n - 1 ] ) } = cs( [ 1, abs( $b ) + 1, $xy + $n ] );
			$n++;
		} while $n < $yz;
	} else 
	{
		$n = 0;
		
		do
		{
			$h_map{ cs( [ 0, 0, $n ] ) } = cs( [ 1, 0, $xy - $n - 1 ] );
            $h_map{ cs( [ 1, 0, $xy - $n - 1 ] ) } = cs( [ 0, 0, $n ] );
            $h_map{ cs( [ 0, abs( $a ) + 1, $n ] ) } = cs( [ 1, abs( $b ) + 1, $x - $n - 1 ]  );
            $h_map{ cs( [ 1, abs( $b ) + 1, $x - $n - 1 ] ) } = cs( [ 0, abs( $a ) + 1, $n ]  );
            $n++;
		} while ( $n < $xy );
		
		$n = 0;
		
		do
		{
            $h_map{ cs( [ 1, 0, $x + $l + $n ] ) } = cs( [ 2, 0, $z - $n - 1 ] );
            $h_map{ cs( [ 2, 0, $z - $n - 1 ] ) } = cs( [ 1, 0, $x + $l + $n ] );
            $h_map{ cs( [ 1, abs( $b ) + 1, $x + $l + $n ] ) } = cs( [ 2, abs( $c ) + 1, $z - $n - 1 ] );
            $h_map{ cs( [ 2, abs( $c ) + 1, $z - $n - 1 ] ) } = cs( [ 1, abs( $b ) + 1, $x + $l + $n ] );
            $n++;
		} while $n < $yz;
		
		$n = 0;
		
		do
		{
            $h_map{ cs( [ 1, 0, $x + $n ] ) } = cs( [ 1, 0, $y - 1 - $n ] );
            $h_map{ cs( [ 1, 0, $y - 1 - $n ] ) } = cs( [ 1, 0, $x + $n ] );
            $h_map{ cs( [ 1, abs( $b ) + 1, $x + $n ] ) } = cs( [ 1, abs( $b ) + 1, $y - 1 - $n ] );
            $h_map{ cs( [ 1, abs( $b ) + 1, $y - 1 - $n ] ) } = cs( [ 1, abs( $b ) + 1, $x + $n ] );
            $n++;
		} while $n < $l;
		
	}
	
	return( %h_map )
}

sub det {
  my @m = @_;
  if (@m == 0) {
    die "should never get here";
  }
  if (@m == 1) {
    return $m[0][0];
  }
  if (@m == 2) {
    my $det = det2(@m);
    return $det;
  }
  my $det = {};
  for my $i (0 .. $#m) {
    my %entry = %{$m[0][$i]};
    next if !%entry;  # skip zero entries

    my @minor = minor($i, @m);
    my $det_minor = det(@minor);

    my $poly_str = poly_str(%$det_minor);

    next if !%$det_minor;
    $det_minor = poly_mult_scalar( (-1) ** $i, $det_minor );
    $det_add = poly_mult(\%entry, $det_minor);
    $det = poly_plus($det, $det_add);
  }
  return $det;
}

sub minor {
# the minor of a matrix.
  my($i, @m) = @_;
  my @minor;
  for my $r (1 .. $#m) {
    my @row = @{$m[$r]};
    my @rowi = grep $_ != $i, 0 .. $#row;
    my @rowm = map $m[$r][$_], @rowi;
    push @minor, [ @rowm ];
  }
  return @minor;
}

sub det2 {
# the determinant of a 2x2 matrix.
  my @m = @_;
  my $plus = poly_mult($m[0][0], $m[1][1]);
  my $minus = poly_mult($m[0][1], $m[1][0]);
  my $det = poly_minus($plus, $minus);

  return $det;
}

sub poly_mult_scalar {
# $a * %poly.
  my($a, $p) = @_;
  my %p = %$p;
  for (keys %p) {
    $p{$_} *= $a;
  }
  return { %p };
}

sub poly_minus {
# %poly1 - %poly2.
  my($p1, $p2) = @_;
  my $p2_minus = poly_mult_scalar(-1, $p2);
  my $poly_minus = poly_plus($p1, $p2_minus);
  return $poly_minus;
}

sub poly_plus {
# %poly1 + %poly2.
  my($p1, $p2) = @_;
  if (!%$p1) { return $p2 }
  if (!%$p2) { return $p1 }
    
  my($exp1_min, $exp1_max) = (sort { $a <=> $b } keys %$p1)[0, -1];
  my($exp2_min, $exp2_max) = (sort { $a <=> $b } keys %$p2)[0, -1];
  my $exp_min = $exp1_min < $exp2_min ? $exp1_min : $exp2_min;
  my $exp_max = $exp1_max > $exp2_max ? $exp1_max : $exp2_max;
  
  my %p;
  for my $exp ($exp_min .. $exp_max) {
    my $coeff = ($p1->{$exp} || 0) + ($p2->{$exp} || 0);
    $p{$exp} = $coeff;
  }

  return { %p };
} 

sub poly_mult {
# %poly1 * %poly2.
  my($p1, $p2) = @_;
  my %p;
  for my $exp1 (keys %$p1) {
    for my $exp2 (keys %$p2) {
      my $exp = $exp1 + $exp2;
      my $coeff = $p1->{$exp1} * $p2->{$exp2};
      $p{$exp} += $coeff;
    }
  }
  return { %p };
}

sub poly_prune {
# remove leading and trailing zeros.
  my $p = shift;
  my %p = %$p;
  my %pruned = %p;
  for my $exp (sort { $a <=> $b } keys %pruned) {
    last if $pruned{$exp};
    delete $p{$exp};
  }
  %pruned = %p;
  for my $exp (sort { $b <=> $a } keys %pruned) {
    last if $pruned{$exp};
    delete $p{$exp};
  }
  return { %p };
}

sub poly_center {
  my $p = shift;
  my %p = %$p;

  my $expc = -1 * (keys(%p) - 1) / 2;
  my %pc;
  for my $exp (sort { $a <=> $b } keys %p) {
    $pc{$expc} = $p{$exp};
    $expc++;
  }
  # make the center term positive, if it isn't.
  if ($pc{0} < 0) {
    for my $exp (keys %pc) { $pc{$exp} *= -1 }
  }

  return { %pc };  
}

sub rolfsen_str {
# "[ 1 - 2 + 1"
  my $p = shift;
  my %p = %$p;
  my $rs;
  my $exp = 0;
  {
    my $coeff = $p{$exp};
    $rs .= $coeff < 0 ? ' - ' : ' + ';
    $rs .= abs($coeff);
    $exp++;
    last if !exists $p{$exp};
    redo;
  }
  $rs =~ s/^ \+ //;
  $rs = '[ ' . $rs;
  return $rs;
}

sub print_matrix {
  my @m = @_;

  # find the longest term length.
  my $length = 0;
  for my $row (@m) {
    my @row = @$row;
    for my $c (0 .. $#row) {
      my %poly = %{$row[$c]};
      my $poly_str = poly_str(%poly);
      if (length($poly_str) > $length) { $length = length($poly_str) }
    }
  }
  $length++;  # for a space

  # print it.
  for my $row (@m) {
    my @row = @$row;
    print "[ ";
    for my $c (0 .. $#row) {
      my %poly = %{$row[$c]};
      my $poly_str = poly_str(%poly);
      printf "%-${length}s", $poly_str;
    }
    print " ]\n";
  }
  return;
}

sub jacobian_matrix {
  my @relators = @_;

  my @m;

  # don't use the last one.
  for my $el1 (0 .. $crossings - 2) {
    for my $el2 (0 .. $crossings - 2) {
      my %poly = differential($relators[$el1], $el2);
      my $poly_str = poly_str(%poly);
      $m[$el1][$el2] = { %poly };
    }
  }
  return @m;
}

sub differential {
# d($ret)/d($el)
  my($ret, $el) = @_;
  my %poly; 
  for my $j (1 .. length($ret)) {
    my $substr = substr($ret, 0, $j);
    next unless my($last) = $substr =~ /(\($el'?\))$/;
    my $coeff;
    if ($last =~ /'/) {
      $coeff = -1;
      $substr =~ s/\Q$last\E$//;
    }
    else {
      $coeff = 1;
    }
    my $plus = () = $substr =~ /(\(\d+\))/g;
    my $minus= () = $substr =~ /(\(\d+'\))/g;
    my $exp = $plus - $minus;
    $poly{$exp} += $coeff;
  }
  return %poly;
}

sub poly_str {
  my %poly = @_;
  return '0' if !%poly;
  my @terms;
  for my $exp (sort { $a <=> $b } keys %poly) {
    my $coeff = $poly{$exp};
    my $term = term_str($coeff, $exp);
    push @terms, $term;
  }
  my $str = join ' ', @terms;
  $str =~ s/^\+ //;
  $str =~ s/^\- /-/;
  return $str;
}

sub term_str {
  my($coeff, $exp) = @_;
  # usual forms.
  my $coeff_str = $coeff >= 0 ? ('+ ' . $coeff) : ('- ' . abs($coeff));
  my $exp_str = 'x^' . $exp;
  if ($exp == 0) {
    $exp_str = '';
  }
  elsif ($exp == 1) {
    $exp_str = 'x';
  }
  else {
  }

  my $str = "$coeff_str$exp_str";
  # " 1x" -> " x"
  $str =~ s/ 1x/ x/;

  return $str;
}

sub relators {
  my %crossings = @_;
  for my $hrc (sort by_hrc keys %crossings) {
    my($el_under, $dir_under) = @{ $crossings{$hrc}{under} };
    my($el_over, $dir_over) = @{ $crossings{$hrc}{over} };

    my $el_next = ($el_under + 1) % $crossings;
    my($exp1, $exp2, $exp_str);
    if (($dir_under eq 'E' && $dir_over =~ /^S/)
     || ($dir_under eq 'W' && $dir_over =~ /^N/)) {
      $exps = $ELS[$el_over];
      $exp1 = "$el_over'";
      $exp2 = $el_over;
    }
    else {
      $exps = uc $ELS[$el_over];
      $exp1 = $el_over;
      $exp2 = "$el_over'";
    }

    # looks better inverted.
    my $relator = "($exp1)($el_next')($exp2)($el_under)";
    push @relators, $relator;

    if ( $DEBUG ) { print "($el_under) = ($el_next)^($exp2)\n"; }
  }
  return @relators;
}

sub print_crossings {
  my %crossings = @_;
  for my $hrc (sort by_hrc keys %crossings) {
    my($el_under, $dir_under) = @{ $crossings{$hrc}{under} };
    my($el_over, $dir_over) = @{ $crossings{$hrc}{over} };
    print "($hrc): under($el_under, $dir_under) over($el_over, $dir_over)\n";
  }
  return;
}


sub path {
  my($t, $i, $h_map) = @_;
  my @t = @$t;
  my @i = @$i;
  my %h_map = %$h_map;

  my @path;

  # the crossings.
  my %crossings;

  # starting stuff.
  my @els;
  my $el = 0;
  push @els, $el;
  my($h, $r, $c) = (0, 0, 0);
  push @path, [$h, $r, $c];
  my $next = $t[$h] == 0 ? 'S' : $t[$h] > 0 ? 'SW' : 'SE';
  if( $DEBUG )
  {
	  print "EL: $el\n";
	  print "PATH: ", cs(@path), " | NEXT($next)\n";
  }

  my $count = 0;
  my $tops = 0;
  PATH: {
    if ($next eq 'N') {
      ($h, $r, $c) = ($h, $r - 1, $c);
    }
    elsif ($next eq 'S') {
      ($h, $r, $c) = ($h, $r + 1, $c);
    }
    elsif ($next eq 'E') {
      ($h, $r, $c) = ($h, $r, $c + 1);
      # hit the right side: go either NW or SW based on t for this handle.
      if ($c == $i[$h]) {
        $next = $t[$h] < 0 ? 'NW' : 'SW';
      }
      # new element unless hit the right side.
      else {
	$crossings{cs($h,$r,$c)}{under} = [ $el, 'E' ];
        if( $DEBUG ) { print "EL: $el -> "; }
        $el++;
	$el = $el % $crossings;
        push @els, $el;
        if( $DEBUG ) { print "$el\n"; }
      }
    }
    elsif ($next eq 'W') {
      ($h, $r, $c) = ($h, $r, $c - 1);
      # hit the left side: go either NE or SE based on t for this handle.
      if ($c == 0) {
        $next = $t[$h] < 0 ? 'SE' : 'NE';
      }
      # new element and a crossing unless hit the left side.
      else {
	$crossings{cs($h,$r,$c)}{under} = [ $el, 'W' ];
        if ( $DEBUG ) { print "EL: $el -> "; }
        $el++;
	$el = $el % $crossings;
        push @els, $el;
        if( $DEBUG ) { print "$el\n"; }
      }
    }
    elsif ($next eq 'SE') {
      # if ending at the bottom, and t < 0 for this handle, need to alter.
      if ($r == abs($t[$h]) && $t[$h] < 0) {
        ($h, $r, $c) = ($h, $r + 1, $c);
      }
      else {
        ($h, $r, $c) = ($h, $r + 1, $c + 1);
      }
      # hit the right side: go W.
      if ($c == $i[$h]) {
        $next = 'W';
      }
      else {
        # a crossing, if not at bottom.
	if ($r <= abs($t[$h])) {
	  $crossings{cs($h,$r,$c)}{over} = [ $el, 'SE' ];
	}
      }
    }
    elsif ($next eq 'SW') {
      # if starting from the top, and t > 0 for this handle, need to alter.
      if ($r == 0 && $t[$h] > 0) {
        ($h, $r, $c) = ($h, $r + 1, $c);
      }
      else {
        ($h, $r, $c) = ($h, $r + 1, $c - 1);
      }
      # hit the left side: go E.
      if ($c == 0) {
        $next = 'E';
      }
      else {
        # a crossing, if not at bottom.
	if ($r <= abs($t[$h])) {
	  $crossings{cs($h,$r,$c)}{over} = [ $el, 'SW' ];
	}
      }
    }
    elsif ($next eq 'NE') {
      # if ending at the top, and t > 0 for this handle, need to alter.
      if ($r == 1 && $t[$h] > 0) {
        ($h, $r, $c) = ($h, $r - 1, $c);
      }
      else {
        ($h, $r, $c) = ($h, $r - 1, $c + 1);
      }
      # hit the right side: go W.
      if ($c == $i[$h]) {
        $next = 'W';
      }
      else {
        # a crossing, if not at top.
	if ($r > 0) {
	  $crossings{cs($h,$r,$c)}{over} = [ $el, 'NE' ];
	}
      }
    }
    elsif ($next eq 'NW') {
      # if starting from the bottom, and t < 0 for this handle, need to alter.
      if ($r > abs($t[$h])) {
        ($h, $r, $c) = ($h, $r - 1, $c);
      }
      else {
        ($h, $r, $c) = ($h, $r - 1, $c - 1);
      }
      # hit the left side: go E.
      if ($c == 0) {
        $next = 'E';
      }
      else {
        # a crossing, if not at top.
	if ($r > 0) {
	  $crossings{cs($h,$r,$c)}{over} = [ $el, 'NW' ];
	}
      }
    }
    else {
      die "unknown NEXT($next)";
    }

    push @path, [$h, $r, $c];
    if( $DEBUG ) { print "PATH: ", cs($h,$r,$c), " | NEXT($next)\n"; }

    # hit the top?
    if ($r == 0) {
      my $h_map_next = $h_map{cs($h, $r, $c)}
        or die "can't find h_map for ($h, $r, $c)";
      ($h, $r, $c) = sc($h_map_next);
      push @path, [$h, $r, $c];
      $next = $t[$h] == 0 ? 'S' : $t[$h] > 0 ? 'SW' : 'SE';
      if( $DEBUG ) { print "PATH: ", cs($h,$r,$c), " | NEXT($next)\n"; }
      $tops += 2;
      if( $DEBUG ) { print "TOPS($tops)\n"; }
    }

    # hit the bottom?
    if ($r > abs($t[$h])) {
      my $h_map_next = $h_map{cs($h, $r, $c)}
        or die "can't find h_map for ($h, $r, $c)";
      ($h, $r, $c) = sc($h_map_next);
      push @path, [$h, $r, $c];
      $next = $t[$h] == 0 ? 'N' : $t[$h] > 0 ? 'NE' : 'NW';
      if( $DEBUG ) { print "PATH: ", cs($h,$r,$c), " | NEXT($next)\n"; }
    }

    last PATH if $h == 0 && $r == 0 && $c == 0;
    $count++; last if $count > 1_000_000;
    redo PATH;
  }

  return(\@els, \@path, $tops, %crossings);
}

sub h_map_str {
# returns a handle map as a single string, for testing.
# "1 - 6, 2 - 3, 4 - 5";
  my %h_map = @_;
  my %c2i;
  my $i = 0;
  for my $from (sort by_hrc keys %h_map) {
    $i++;
    $c2i{$from} = $i;
  }
  my @str;
  my %seen;
  for my $from (sort by_hrc keys %h_map) {
    next if $seen{$from};
    my $to = $h_map{$from};
    $seen{$to}++;
    my $from_i = $c2i{$from};
    my $to_i = $c2i{$to};
    push @str, "$from_i - $to_i";
  }
  my $str = join ', ', @str;
  return $str;
}

sub by_hrc {
  my @hrca = sc($a);
  my @hrcb = sc($b);
  $hrca[0] <=> $hrcb[0]
    or
  $hrca[1] <=> $hrcb[1]
    or
  $hrca[2] <=> $hrcb[2]
    ;
}

sub print_map {
  my %map = @_;
  for my $from (sort by_hrc keys %map) {
    print "$from <-> $map{$from}\n";
  }
  return;
}

sub cs {
# the coordinate-string for ($h, $r, $c).
  my $hrc = shift;
  my @hrc;
  if (ref($hrc) eq 'ARRAY') {
    @hrc = @$hrc;
  }
  else {
    @hrc = ($hrc, @_);
  }
  local $" = ', ';
  my $cs = "@hrc";
  return $cs;
}

sub sc {
# coord-string -> ($h, $r, $c)
  my $hrc = shift;
  my @sc = split /,\s*/, $hrc;
  return @sc;
}

sub verify_i {
  my @i = @_;
  
  # each must be non-negative integers.
  # each must be positive integers.
  my $sum = 0;
  for my $i (@i) {
    ($i >= 0) && (int($i) == $i)
      or error( "each i (@i) must be non-negative." );
    $sum += $i;
  }
  
  # the sum must be even.
  $sum % 2 == 0 or error( "sum of i (@i) must be even." );

  # the middle must be the maximum.
  ($i[1] >= $i[0]) && ($i[1] >= $i[2])
    or error( "middle i (@i) must be the maximum." );

  if ($i[1] > ($i[0] + $i[2])) {
    print "RS: loops\n";
  }

  # all ok.
  return;
}

sub error {
# XXX check this.
# get here if $els->[-1] != 0.  (i.e. did not close the path)
  my($message, $els, $path, %crossings) = @_;
  print $message;
  #print "RS: unknot, link, or unlink\n";
  exit;
}