#!/usr/local/bin/perl -w
# GHG.p                      see also GH.p written april 1998.
#                           this is GHG.p, written 0898, which 
#                           uses Gallager's method of construction.
#
# - which avoids having any cycles of length 4.
# except in version 1, if Ps are superposed, it doesn't do this!
# version 2 980906 fixes this.
# this version (2) used until july 2000.
# Move to version 3 made when investigating staircase codes.
# worked all day on it but couldn't make the stairballoon do its job
#
# Definitely works ok for non-superposed permutation matrices, but I 
# can't remember whether it does the right thing when they are superposed.
#
# used to make GHmaker.p which makes alist matrices based on permutations.
# now does it itself
#
# bugs: it should check to see if you are asking for out-of-range
#       locations
#
# other similar programs:  Aaugment.p Acleanup.p Ashorten.p Ashorten2.p
#                           productA.p
#
# see ~mackay/pub/gallager/GH.html or http://www.inference.phy.cam.ac.uk/mackay/gallager/GH.html 
#
# modified 98 11 02 so that it gives up immediately if failure
#   to make a perm.
# repeats from start if tries > 1
#
# usage:
# this one works:
# GHG.p greedyfirst=1 N=16 verbose=1 spec1 > ../GHC/tmp
# checking version 2: 980906
# GHG.p greedyfirst=1 N=10 verbose=1 spec2 > ../GHC/tmp
# wow, this seems to work! a 10 by 10 superposition of three matrices!
# these don't
#  GHG.p N=24 verbose=1 spec3 > ../GHC/tmp
#  GHG.p N=20 verbose=1 spec1 > ../GHC/tmp
#  GHG.p N=198  spec3 > ../GHC/tmp
#  GHG.p N=999  spec9 > ../GHC/tmp
# almost worked ^^
#  GHG.p N=999 seed=234 spec9 > ../GHC/tmp2
# 4000 bits is permitted in disk drives
#  GHG.p N=3996 seed=1 spec9 > ../GHC/seagate4
# this worked!
#  GHG.p N=1998 seed=1 spec9 > ../GHC/seagate2
# this worked too.
#  GHG.p N=999 seed=1 spec9 > ../GHC/seagate1
# this gets stuck.
#  GHG.p N=999 seed=1 spec99 > ../GHC/seagate1
#  GHG.p N=100 seed=1 specstair2 > ../GHC/stair2
#this works!
#
# see also GHdraw.p
#
# obsolete:
$unthicken=1; # whether to unthicken again, when writing, if thickness > 2
$stairthickness=2 ; # this (2) gives standard staircase.
#
$stairballoon = 0 ; # whether, when making staircases, to
# force all permutation matrices to look out for not only
# 4-cycles, but near 4-cycles.

$forbid_four_cycles = 1 ; # the standard setting.

$writeanyway = 0 ; # whether to write a half-baked imperfect code.
$tries = 1 ;      # number of attempts to make the whole code. 
$greedyfirst = 0 ; # whether to try greedy permutation instead of random
$nattempts = 100 ; # number of timesto try to make a perm before giving up.
$seed="123";
$verbose = 0 ; 
$N=10000 ;
# $argvcommand = join(' ' , @ARGV) ;

eval "\$$1=\$2" while @ARGV && $ARGV[0]=~ /^(\w+)=(.*)/ && shift;

$Norig = $N ; 
print STDERR "===== GHG.p ===== (c) DJCM april 1998 =====\n" ;
$occn = 0 ; $occm = 0 ; # ugly global counters.
#
# read the spec file
#
$status = 0 ; $n = 0 ; $m = 0 ; $lambda = 0 ;
$c = 0 ; # command number
$dlist = ":" ; $dprod = 1 ; 
while(<>) {
    s/^\s*// ; 
    if (/^#/) {
#	print STDERR ;
	next ;
    }
    s/\#.*// ;		#
    if (/lambda/) {
	print STDERR " lambda = " ;
	($n,$m) = /lambda\s+(.*)\/(.*)/ ;
	print STDERR "$n / $m \n" ;
	if ( $m <= 0 ) { $status -- ; }
	else {
	    $lambda = ( 1.0 * $n ) / ( 1.0 * $m ) ; 
	}
	&denoms_add ( $m ) ; 
    } elsif  (/refresh/) {
	$c ++ ;
	$type[$c] = "refresh" ; 
	print STDERR "refresh\n" ;
    } elsif  (/forbidI/) {
	$c ++ ;
	$type[$c] = "forbidI" ;
	($sizen[$c],$sized[$c])  =  /forbidI\s+(.*)\/(.*)/ ;
	print STDERR "forbid permutations like I $sizen[$c]/$sized[$c]\n" ;
	&denoms_add ( $sized[$c] ) ; 
    } elsif  ( /[SIP].*at/ ) {
	$c ++ ;
	($cols[$c],$rows[$c],$type[$c],$sizen[$c],$sized[$c],
	 $xn[$c],$xd[$c],$yn[$c],$yd[$c]) =
	     /(.*)x(.*)\s+([SIP])\s+(.*)\/(.*)\s+at\s+(.*)\/(.*),(.*)\/(.*)/ ;
	print STDERR " $cols[$c] x $rows[$c] $type[$c] size $sizen[$c]/$sized[$c] at $xn[$c]/$xd[$c],$yn[$c]/$yd[$c]\n" ;
	&denoms_add ( $sized[$c] ) ; 
	&denoms_add ( $xd[$c] ) ; 
	&denoms_add ( $yd[$c] ) ;
    }
}
sub denoms_add { # should really check for divisibility here!
    local ($d) = @_ ;
    if (!($dlist =~ /:$d:/) ) {
	$dlist = $dlist."$d:" ;
	$dprod = $dprod * $d ;
    }
}
$C = $c ;
print STDERR "# Read $C commands\n" ;
print STDERR "# suggest Make M divisible by $dprod\n" ;
	
# suggested value for $M
$M = $N * $m / $n ;
$M = ( int ( $M / $dprod ) ) * $dprod ;
$N = $M * $n / $m ;
print STDERR "# suggest M = $M, N = $N\n" ;
if ( $N != $Norig ) {
    $N = $Norig ;
    $M = int ( $N * $m / $n ) ;
    $N = $M * $n / $m ;
    print STDERR "WARNING WARNING!!!!!!!!!!!!!!!!!!!!\n" ; 
    print STDERR "# now running with \n" ;
    print STDERR "# M = $M, N = $N\n" ;
}
print STDERR "\nP1\n$N $M\n" ; # useful for gv 

$codeattempt = 0 ; 
# repeat until successful or run out of tries.
MAINLOOP: do {
    $stuck = 0 ; 
    $codeattempt ++ ;
    if ( $codeattempt > 1 ) { $seed ++ ; print STDERR "seed=$seed\n" ; }
    srand ( $seed ) ;
    
    $head = $N." ".$M."\n" ;

    undef(@mntaken) ; undef(@mnforbidden) ; # clear, delete, forget, cancel
    undef(%mntaken) ; undef(%mnforbidden) ; # clear, delete, forget, cancel
# make null lists 
    for ( $n = 1 ; $n <= $N ; $n ++ ) {
	$nlist[$n] = "" ; 
	$numnlist[$n] = 0 ;
    }
    for ( $m = 1 ; $m <= $M ; $m ++ ) {
	$mlist[$m] = "" ; 
	$nummlist[$m] = 0 ;
    }
    
    $failed=0;
    &refresh_perm() ;

    for ( $c = 1 ; (!($failed)) && $c <= $C ; $c ++ ) {
	$thetype = $type[$c] ;
	if ( $type[$c] eq "refresh" ) {
	    &refresh_perm() ; 
	} elsif ( $type[$c] eq "I" ) {
	    &do_I($cols[$c],$rows[$c],$sizen[$c],$sized[$c],
		  $xn[$c],$xd[$c],$yn[$c],$yd[$c]) ;  
	} elsif ( $type[$c] eq "S" ) {
	    &do_S($cols[$c],$rows[$c],$sizen[$c],$sized[$c],
		  $xn[$c],$xd[$c],$yn[$c],$yd[$c]) ;  
	} elsif ( $type[$c] eq "forbidI" ) {
	    &registerI($M*$sizen[$c]/$sized[$c]) ; 
	} elsif ( $type[$c] eq "P" ) {
	    &do_P($cols[$c],$rows[$c],$sizen[$c],$sized[$c],
		  $xn[$c],$xd[$c],$yn[$c],$yd[$c]) ;
	} else {
	    print STDERR "oops, unknown type $type[$c]\n"  ;
	}
    }
    print STDERR "\n" ; 
    print STDERR "\n failures: $failed.\n" ; 

} while ( ($failed) && ( $codeattempt < $tries ) ) ;


if ( $writeanyway || ($failed==0) ) {
# write code

    # first remove any thickenning from the stair case
    if ( $unthicken && ($stairthickness > 2 ) ) {
	for ( $c = 1 ; $c <= $C ; $c ++ ) {
	    if ( $type[$c] eq "S" ) {
		&undo_S($cols[$c],$rows[$c],$sizen[$c],$sized[$c],
			$xn[$c],$xd[$c],$yn[$c],$yd[$c]) ;
	    }
	}
    }
    
    print STDERR "Writing code.\n" ;

# find $t (max) and $tr (max)

    $t = 0 ; $tr = 0 ; 
    for ( $n = 1 ; $n <= $N ; $n ++ ) {
	if ( $numnlist[$n] > $t ) {
	    $t = $numnlist[$n] ;
	}
    }
    for ( $m = 1 ; $m <= $M ; $m ++ ) {
	if ( $nummlist[$m] > $tr ) {
	    $tr = $nummlist[$m] ;
	}
    }
    $head .= "$t $tr\n" ;
# append zeros to short lines
    for ( $n = 1 ; $n <= $N ; $n ++ ) {
	$thisn = $numnlist[$n] ;
	$head .= "$thisn " ; 
	while ( $thisn < $t ) {
	    $nlist[$n] .= "0\t" ;
	    $thisn ++ ;
	}
    }
    $head .= "\n" ;
    
    for ( $m = 1 ; $m <= $M ; $m ++ ) { 
	$thism = $nummlist[$m] ; 
	$head .= "$thism " ; 
	while ( $thism < $tr ) {
	    $mlist[$m] .= "0\t" ;
	    $thism ++ ;
	}
    }
    $head .= "\n" ;
    
    print $head ; 
# finish lists 
    for ( $n = 1 ; $n <= $N ; $n ++ ) {
	$nlist[$n] =~ s/\t$/\n/ ; 
	print $nlist[$n] ; 
    }
    for ( $m = 1 ; $m <= $M ; $m ++ ) {
	$mlist[$m] =~ s/\t$/\n/ ; 
	print $mlist[$m] ; 
    }
}
# print STDERR  "forbidden = $#mnforbidden\n"  ; 

# define some useful subroutines
sub refresh_perm {
    $nrandperms = 0 ;
}
# make a staircase
sub do_S {
    local($Co,$Ro,$sn,$sd,$xn,$xd,$yn,$yd) = @_ ;
    local($n0) = 1 + $M * $xn / $xd ; 
    local($m0) = 1 + $M * $yn / $yd ;
    local($w) = $M * $sn / $sd ; 
    for ( $co = 1 ; $co <= $Co ; $co ++ ) { # run through squares by col
	local($m) = $m0 ; 
	for ( $ro = 1 ; $ro <= $Ro ; $ro ++ ) { # # run through squares by row
	    local($n) = $n0 ; 
	    # now process each square
	    for ( $ww = 1 ; $ww <= $w ; $ww ++ , $n ++ , $m ++ ) {
		&mnadd () ;
	    }
	}
	local($m) = $m0 ; # back to the top, now do all rows except last
	for ( $ro = 1 ; $ro <= $Ro ; $ro ++ ) {# run through squares by row
	    local($n) = $n0 + $w - 1 ; # top right hand corner wrap-around
	    &mnadd () ; $m ++ ; 
	    local($n) = $n0 ; 
	    for ( $ww = 2 ; $ww <= $w ; $ww ++ , $n ++ , $m ++ ) {
		&mnadd () ;
	    }
	}
	if ($stairthickness > 2) { # add some more 1s to effectively
				# forbid longer cycles that use the trellis
	    print STDERR "Thickenning staircase\n" ;
	    $thicken = 1 ; 
	}
	for ($sss = 3  ; $sss <= $stairthickness ; $sss += 2 ) {
	    $thicken ++; # 2, on first time round
	    local($m) = $m0 ; # back to the top
	    for ( $ro = 1 ; $ro <= $Ro ; $ro ++ ) {
		local($n) = $n0 + $w - $thicken ;
		for ( ; $n < $n0 + $w ; $n ++ , $m ++ ) {
		    &mnadd () ; 
		}
		local($n) = $n0 ; 
		for ( $ww = $thicken+1 ; $ww <= $w ; $ww ++ , $n ++ , $m ++ ) {
		    &mnadd () ;
		}
	    }			# 
	    #    and now thicken the far side (upper) too.
	    local($m) = $m0 ; # back to the top
	    for ( $ro = 1 ; $ro <= $Ro ; $ro ++ ) {
		local($n) = $n0 + $thicken - 1 ;
		for ( ; $n < $n0 + $w ; $n ++ , $m ++ ) {
		    &mnadd () ; 
		}
		local($n) = $n0 ; 
		for ( $ww = -$thicken+2 ; $ww <= 0 ; $ww ++ , $n ++ , $m ++ ) {
		    &mnadd () ;
		}
	    }			# 
	}			# 
	$n0 += $w;		# 
    }				# 
    if ($verbose) { 
	print STDERR "\n" ; 
    }
}
sub do_I {
    local($Co,$Ro,$sn,$sd,$xn,$xd,$yn,$yd) = @_ ;
    local($n0) = 1 + $M * $xn / $xd ; 
    local($m0) = 1 + $M * $yn / $yd ;
    local($w) = $M * $sn / $sd ; 
    for ( $co = 1 ; $co <= $Co ; $co ++ ) {
	local($m) = $m0 ; 
	for ( $ro = 1 ; $ro <= $Ro ; $ro ++ ) {
	    local($n) = $n0 ; 

	    for ( $ww = 1 ; $ww <= $w ; $ww ++ , $n ++ , $m ++ ) {
		&mnadd () ;
	    }
	}
	$n0 += $w;
    }				# 
    # each time I is called,
    # register the identity as a randperm
    if ($verbose) { 
	print STDERR "\n" ; 
    }
}


# Gallager's method for P. (replaces old_do_P)
sub do_P {
    local($Co,$Ro,$sn,$sd,$xn,$xd,$yn,$yd) = @_ ;
    local($n0) =  $M * $xn / $xd ; 
    local($m0) =  $M * $yn / $yd ; 
    local($w) = $M * $sn / $sd ; 
    local($n) = $n0 ; 
    $valid = 1 ; 
    for ( $co = 1 ; $valid && $co <= $Co ; $co ++ , $n += $w ) {
	local($m) = $m0 ; 
	for ( $ro = 1 ; $valid && $ro <= $Ro ; $ro ++ , $m += $w ) {

	    &invent_a_perm() ;
	    if (!$valid) { $failed ++ ; } 
	}
    }
    
}

sub create_eligibility_square {
    for ( $nq = 1 ; $nq <= $w ; $nq ++ ) {
	$numeligiblem[$nq] = 0 ;
	$numeligiblen[$nq] = 0 ;
    }
    for ( $nq = 1 ; $nq <= $w ; $nq ++ ) {
	  local($nn) = $n + $nq ;
	  for ( $mq = 1 ; $mq <= $w ; $mq ++ ) {
	      local($mm) = $m + $mq ;
	      $mn = "$mm,$nn" ; # true m n coordinates
	      $mnq = "$mq,$nq" ; # local coordinates
	      if ( $forbid_four_cycles && ( ( $mnforbidden{$mn} || $mntaken{$mn}) ) ) {
		  $eligible{$mnq} = 0 ;
		  if ( $mntaken{$mn}) {
		      $eligible{$mnq} = -3 ; # this means already taken
		  }
	      } else {
		  $eligible{$mnq} = 1 ;
		  $numeligiblen[$nq] ++ ;
		  $numeligiblem[$mq] ++ ;
	      }
	  }
      }
    if ( $verbose >=2 || ($verbose&& $a==1) ) { &show_square () ; }
    for ( $q = 1 ; $q <= $w ; $q ++ ) {
	if ( !$numeligiblem[$q] ||   	!$numeligiblen[$q] ) {
	    print STDERR "definitely stuck! $mn $q " ;
	    print STDERR "$numeligiblem[$q] $numeligiblen[$q] \n" ;
#	    exit(0) ;
	    $stuck = 1 ;  $failed ++ ;
#	    next MAINLOOP ;
	}
    }
}
sub invent_a_perm{
    if ( $verbose > 0 ) {
	printf STDERR "\n== Starting perm %d %d == " , $nrandperms+1 , $w ; 
    } else {
	printf STDERR " P%d %d " , $nrandperms+1 , $w ; 
    } 

    # here we use a two-pass procedure. 
# in the first pass, we make a tentative addition of each m,n
# in the second we go through adding them to the master list.
    
  ALOOP: for ( $a  = 1 ; $a <= $nattempts ; $a ++ )  { 
      print STDERR "p" ;
      &create_eligibility_square() ;
#  Run through columns, picking rows.
      if ( !$stuck ) {

	  $valid = 1 ; # nq is the column number in relative coordinates
	QLOOP:
	  for ( $nq = 1 ; $nq <= $w ; $nq ++ ) {
	      if (  $numeligiblen[$nq] >= 1 ) {
		  if ( $a > 1 || !$greedyfirst )  { # pick a random eligible row.
		      $q = int ( rand( $numeligiblen[$nq] ) ) + 1 ;
		  } else { # use greedy choice.
		      $q = 1 ; 
		  }
	      } else {
		  if ( $verbose > 0 ) {
		      print STDERR "stuck row $nq of $w\n" ;
		  } else {
#		  print STDERR "s" ; 
		  }
		  $valid = 0 ; 
		  last QLOOP ; # at this point, could invoke the Gallager emergency 
		  # procedure.
	      }
	      if ( $verbose>=2) {
		  print STDERR "($q/$numeligiblen[$nq])" ;
	      }
	      local($counter) =  0 ; 
	      # we are still within QLOOP. Here we clumsily identify 
	      # the q'th eligible row.
	      
	      $p[$nq] = 0 ;
	    FLOOP:		
	      for ( $mq = 1 ; $mq <= $w ; $mq ++ ) {
		  $mnq = "$mq,$nq" ;
		  if ($eligible{$mnq}>0) {
		      $counter ++ ;
		  }
		  if ($counter == $q ) {
		      $p[$nq] = $mq ; # permutation entry
		      last FLOOP ;
		  }
	      } # end of FLOOP
# we have now made a tentative valid selection for this row of permutation
	      if ( $p[$nq] ) { # update the eligibility square, but not yet the 
# master lists.
		  if ($verbose) { 
		      print STDERR " $q) $nq>$mq; " ;
		  }
# occn and occm are reserved global variables for counting occupancy
		  if ( $occn > 0 ) {
		      $occn = 0 ;
		  }
		  if ( $occm > 0 ) {
		      $occm = 0 ;
		  }
		  for ( $nq2 = 1 ; $nq2 <= $w ; $nq2 ++ ) {
		      $mnq2 = "$mq,$nq2" ; #  (*) 
		      if ( $nq2 == $nq ) {
			  $eligible{$mnq2} += 10 ; # mark the new '1'
		      } elsif ( $eligible{$mnq2} > 0 ) {
			  $numeligiblen[$nq2] -- ;
			  $eligible{$mnq2} *= -1 ; # remove 1s in this row/col
		      } elsif ( $eligible{$mnq2} == -3 ) {
			  $alreadyoccn[$occn] = $nq2 ; $occn ++ ;
		      }
		  }
		  # now do everyone in our column
		  for ( $mq2 = 1 ; $mq2 <= $w ; $mq2 ++ ) {
		      $mnq2 = "$mq2,$nq" ;
		      if ( $eligible{$mnq2} == -3 ) {
			  $alreadyoccm[$occm] = $mq2 ; $occm ++ ;
			  if ( $stairballoon ) { # also warn about neighbours
			      for ( $iii = -$stairballoon ; $iii <= $stairballoon ; $iii ++ ) {
				  if ($iii == 0 ) { $iii ++ ; }
				  $mq22 = $mq2 + $iii ; 
				  if ($mq22 < 1 ) { $mq22 += $w ; }
				  if ($mq22 > $w ){ $mq22 -= $w ; }
			      }
			      $alreadyoccm[$occm] = $mq22 ; $occm ++ ;
			  }
		      }
		  }
		  if ( $stairballoon ) { # go to the whole rest of matrix too.
		      local($nn) = $n + $nq ;
		      local(@m2list) = split (  /\s+/ , $nlist[$nn] ) ;
		      # this gives us ^ all the rows that this 1 is
		      # connected to. We are intersted in all these rows +/- several
		      foreach $m2 ( @m2list ) {
			  for ( $iii = -$stairballoon ; $iii <= $stairballoon ; $iii ++ ) {
			      $mq2 = $m2 + $iii ;    
			      if ($mq2 < 1 ) { $mq2 += $M ; }
			      if ($mq2 > $M ){ $mq2 -= $M ; }
			      # go and get the implicated columns, which we are worried about
			      local(@n2list) = split ( /\s+/ , $mlist[$mq2] ) ;
			      foreach $n2 ( @n2list ) {
				  $nnnn = $n2 - $n ;
				  if (($nnnn <= $w) && ($nnnn>=1)) {
				      $alreadyoccn[$occn] = $nnnn ; $occn ++ ;
#				      print STDERR "phew, banned something!\n" ; 
				  }
			      }
			  }
		      }
		  }
		  if ( $stairballoon ) { # also set $mq to other nearby values
		      # to establish whether nearby rows are occupied
IIILOOP:	      for ( $iii = -$stairballoon ; $iii <= $stairballoon ; $iii ++ ) {
			  $mq2 = $mq + $iii ; 				# 
			  if ($mq2 < 1 ) { $mq2 += $w ; }
			  if ($mq2 > $w ){ $mq2 -= $w ; }
			  $alreadyoccm[$occm] = $mq2 ; $occm ++ ;
			  if ($iii == 0 ) { next IIILOOP ; }
			  # OK, let's do the above stuff (*) again.
			  for ( $nq2 = 1 ; $nq2 <= $w ; $nq2 ++ ) {
			      $mnq2 = "$mq2,$nq2" ;
			      if ( $eligible{$mnq2} == -3 ) {
				  $alreadyoccn[$occn] = $nq2 ; $occn ++ ;
			      }
			  }
		      }
		  }
# version 2: (for cases where matrices are superposed)
# Need also to remove from 		  $eligible{$mnq} 
# any points that could now make a square. - as in mnadd
#
# Here it is sufficient to identify all 1s in this col and this row and outer
# product them.
		  if (( $occn > 0 )&&( $occm > 0 )) {
		      $nnnn = $occn * $occm ; 
		      print STDERR "inel'ing $nnnn ; " if ($verbose>4);
		      foreach $mq2 (@alreadyoccm[0..$occm-1]) {
			  foreach $nq2 (@alreadyoccn[0..$occn-1]) {
			      $mnq2 = "$mq2,$nq2" ;
			      print STDERR "$mnq2 " if($verbose>2); 
#			      if ( -e $eligible{$mnq2} ) {
			      if ( $eligible{$mnq2} > 0 ) {
				  $numeligiblen[$nq2] -- ;
				  $eligible{$mnq2} *= -1 ;
			      }
#			      } else { print STDERR "$mnq2 not exist\n" ; }
			  }
		      }
		  }
# end ineligibilification . 
		  if ( $verbose >= 2 ) { # show square
		      print STDERR "New square:" ;
		      &show_square () ;
		  }
	      }  else  {
		  if ($verbose) { 
		      print STDERR " $p[$nq] - permutation generator stuck row $nq of $w\n" ;
		  }
		  $valid = 0 ; 
		  last QLOOP ;
	      }
	  } # end of QLOOP
      } else { # if we are stuck
	  $valid = 0 ; 
	  last ALOOP ;
      }
      if ( !$valid ) {
	  if ($verbose>=2) { print STDERR "start again\n" ; }
      } else {
	  last ALOOP ;
      }
  }
    if ( !$valid ) {
	if ($verbose>=0) { print STDERR "\nfailed to make perm $w,$n,$m\n" ; }
    } else {
	&permadd($w,$n,$m);
print STDERR "a" ;				# 
	if ( $verbose >=1 ) {
	    &show_all () ;
	}
	$nrandperms ++ ;
    }
}

sub show_square {
    print STDERR "\n" ;
    for ( local($mq) = 1 ; $mq <= $w ; $mq ++ ) {
	for (local( $nq) = 1 ; $nq <= $w ; $nq ++ ) {
	    local($mnq) = "$mq,$nq" ;
	    printf STDERR "%2d " , $eligible{$mnq} ;
	}
	print STDERR "\n" ;
    }    
}

sub show_all {
    print STDERR "\n" ;
    for ( local($mq) = 1 ; $mq <= $M ; $mq ++ ) {
	for (local( $nq) = 1 ; $nq <= $N ; $nq ++ ) {
	    local($mn) = "$mq,$nq" ;
	    $nosir = " " ;
	    if ( $mnforbidden{$mn} ) { $nosir = "." ; }
	    if ( $mntaken{$mn} ) { $nosir = 1 ; }
	    printf STDERR "%1s " , $nosir ; 
	}
	print STDERR "\n" ;
    }    
}

sub mnadd {
    local($mn)="$m,$n";
    if ($mnforbidden{$mn}) { print STDERR "warning, $mn forbidden!\n" ; }
    if ($verbose) { 
	print STDERR "add $mn " ; 
	if ( $verbose >=2 ) {
	    print STDERR "\n$m: splitting >$mlist[$m]<\n" ; 
	}
    }
    if ( $mntaken{$mn} ) {
	print STDERR "ERR1: already taken m,n = $mn ($M,$N)\n"; 
    }
    $mntaken{$mn} = 1 ;
    # now handle prevention of four-cycles
    # original version only handles four-cycles (in stages (1)-(3))
    # but now (000704) am adding extra capability to kill 6 or eight-cycles
    # for special case of staircase.
    
    # (1)look at all on this row, find their vertical connections, and forbid.
    # also look at all on neighbouring rows and do same.
#    if ( ($thetype eq "P") && $stairballoon ) {
#	print STDERR  "yes! " ; 
#    } else {
#	print STDERR "tt= $thetype : " ;  
#    }
    if ( ($thetype eq "P") && $stairballoon ) {
	$combinedmlists = $mlist[$m] ;
	$mm = $m ; $mmm = $m ; 
	for ($bbb = 1 ; $bbb <= $stairballoon ; $bbb ++ ) {
	    $mm ++  ; if ($mm > $M){ $mm = 1 ;}#assuming full height staircase
	    $mmm --  ; if ($mmm < 1){ $mmm = $M ;}#assuming full height staircase
	    $combinedmlists .= " ".$mlist[$mm]." ".$mlist[$mmm] ;		#
	}
	# fake command, to check something happns!
	local(@n2list) = split ( /\s+/ , $mlist[$m] ) ;
	$tmpo =  $#n2list+1 ;
	# end fake 
	local(@n2list) = split ( /\s+/ , $combinedmlists );
	# more fake stuff......
	$tmpo2 =  $#n2list+1 ;
	if ( $tmpo2 <= $tmpo ) {
	    print STDERR "oi, enlarged list same size? $combinedmlists : $m\n" if ($verbose >-1) ;
	}	else {
#	    print STDERR "$tmpo2 > $tmpo : " ; 
	} # end fake
    } else { # standard method
	local(@n2list) = split ( /\s+/ , $mlist[$m] ) ;
	if (  $#n2list+1 !=  $nummlist[$m] ) {
	    print STDERR "ERR2: gave rise to $#n2list+1 entries" ;
	    print STDERR ", cf nummlist = $nummlist[$m]\n" ; 
	}
    } # we now have a list of n's that the current column m is closely
     # connected to
    foreach $n2 ( @n2list ) { 
	print STDERR " $n2: now splitting >$nlist[$n2]<\n" if ($verbose>=2) ;
	local(@m2list) = split (  /\s+/ , $nlist[$n2] ) ;
	if (  $#m2list+1 !=  $numnlist[$n2] ) {
	    print STDERR "ERR3:  gave rise to $#m2list+1 entries, cf numnlist = $numnlist[$n2]\n" ; 
	}

	foreach $m2 ( @m2list ) { 
	    if ( ($thetype eq "P") && $stairballoon ) {
		$mm = $m2 ; $mmm = $m2 ; 
		for ($bbb = 1 ; $bbb <= $stairballoon ; $bbb ++ ) {
		    $mm ++ ; if ($mm > $M) { $mm = 1 ;} #assuming full height staircase
		    $mmm --  ; if ($mmm < 1){ $mmm = $M ;} 
		    $m2n="$mm,$n" ; $mnforbidden{$m2n} = 1 if ($mm  != $m ) ;
		    $m2n="$mmm,$n"; $mnforbidden{$m2n} = 1 if ($mmm != $m ) ;
		}
	    }
	    if ( $m2 != $m ) {
		$m2n="$m2,$n";
		if ( $verbose >=2 ) { print STDERR " forbidding m,n = $m2n ($M,$N)\n"; 	}
		$mnforbidden{$m2n} = 1 ;
	    }
	}
    }
    # (2) look at all on this col, find their horizontals, and forbid loop.
    local(@m2list) = split (  /\s+/ , $nlist[$n] ) ;
    if (  $#m2list+1 !=  $numnlist[$n] ) {
	print STDERR "ERR3a:  gave rise to $#m2list+1 entries, cf numnlist = $numnlist[$n2]\n" ; 
    }
    foreach $m2 ( @m2list ) { 
	local(@n2list) = split ( /\s+/ , $mlist[$m2] ) ;
	if (  $#n2list+1 !=  $nummlist[$m2] ) {
	    print STDERR "ERR2a: gave rise to $#n2list+1 entries, cf nummlist = $nummlist[$m]\n" ; 
	}

	foreach $n2 ( @n2list ) { 
	    if ( $n2 != $n ) {
		$m2n="$m,$n2";
		if ( $verbose >=2 ) {
		    print STDERR " forbidding m,n = $m2n ($M,$N)\n"; 
		}
		$mnforbidden{$m2n} = 1 ;
	    }
	}
    }
    # and (3) look at outer product of my row and my column
    local(@m2list) = split (  /\s+/ , $nlist[$n] ) ;
    local(@n2list) = split ( /\s+/ , $mlist[$m] ) ;
    foreach $m2 ( @m2list ) { 
	foreach $n2 ( @n2list ) { 
	    $m2n="$m2,$n2";
	    if ( $verbose >=2 ) {
		print STDERR " forbidding m,n = $m2n ($M,$N)\n"; 
	    }
	    $mnforbidden{$m2n} = 1 ; 
	}
    }

    # put m and n on each others lists 
    $nlist[$n] .= $m."\t" ;
    $mlist[$m] .= $n."\t" ;
    $numnlist[$n] ++ ;
    $nummlist[$m] ++ ;
    if ( $verbose >=2 ) {
	&show_all () ;
    }
}

sub mnsubtract {
#    local($mn)="$m,$n";
#    local($m,$n)=@_;

    $count = ( $mlist[$m] =~ s/\b$n\t// ) ;
    if ( $count ) {
        $nummlist[$m] -- ;
    } else {
        print STDERR "error, unable to remove $n from mlist $m\n $mlist[$m]\n" ; 
    }
    $count = ( $nlist[$n] =~ s/\b$m\t// ) ;
    if ( $count ) {
        $numnlist[$n] -- ;
    } else {
        print STDERR "error, unable to remove $m from nlist $n\n $nlist[$n]\n" ; 
    }
}

sub permadd {
    local ( $Nl , $noff , $moff ) = @_ ; 
    for ( $nl = 1 ; $nl <= $Nl ; $nl ++ ) {
	local($n) = $noff + $nl ;
	local($m) = $p[$nl] + $moff ; 
	if ( ( $m > $M ) ||  ( $n > $N ) ) {
	    print STDERR "invalid m,n = $m,$n ($M,$N)\n" ;
	} else {
	    &mnadd () ;
	}
    }
}

# unthicken a staircase
sub undo_S {
    local($Co,$Ro,$sn,$sd,$xn,$xd,$yn,$yd) = @_ ;
    local($n0) = 1 + $M * $xn / $xd ; 
    local($m0) = 1 + $M * $yn / $yd ;
    local($w) = $M * $sn / $sd ; 
    for ( $co = 1 ; $co <= $Co ; $co ++ ) { # run through squares by col
	print STDERR "UN-Thickenning staircase\n" ;
	$thicken = 1 ; 
	for ($sss = 3  ; $sss <= $stairthickness ; $sss += 2 ) {
	    $thicken ++; # 2, on first time round
	    local($m) = $m0 ; # back to the top
	    for ( $ro = 1 ; $ro <= $Ro ; $ro ++ ) {
		local($n) = $n0 + $w - $thicken ;
		for ( ; $n < $n0 + $w ; $n ++ , $m ++ ) {
		    &mnsubtract () ; 
		}
		local($n) = $n0 ; 
		for ( $ww = $thicken+1 ; $ww <= $w ; $ww ++ , $n ++ , $m ++ ) {
		    &mnsubtract () ;
		}
	    }			# 
	    #    and now thicken the far side (upper) too.
	    local($m) = $m0 ; # back to the top
	    for ( $ro = 1 ; $ro <= $Ro ; $ro ++ ) {
		local($n) = $n0 + $thicken - 1 ;
		for ( ; $n < $n0 + $w ; $n ++ , $m ++ ) {
		    &mnsubtract () ; 
		}
		local($n) = $n0 ; 
		for ( $ww = -$thicken+2 ; $ww <= 0 ; $ww ++ , $n ++ , $m ++ ) {
		    &mnsubtract () ;
		}
	    }			# 
	}			# 
	$n0 += $w;		# 
    }				# 
    if ($verbose) { 
	print STDERR "\n" ; 
    }
}

##################### END #################