#!/usr/bin/perl
#
# file:         lineq.pl
# purpose:      generates and executes GAP code for solving linear equation on 
#               algebras
# created:      pasha dec 17 2012,mar 22-23,apr 6-8,10-13,15-16,18 2014
# modified:     pasha aug 17-18,nov 11 2020,jul 5 2021
# modification: - put version to 0.4.3
#               - cleanup
#               - ref to the web in --help
#               - when eliminating zeros, allow terms like 0*F(...)
#               - revert broken transformation from 0.4.2
#
# See https://www1.osu.cz/~zusmanovich/soft/lineq/
#
# TODO: 
# - take into account (anti)commutativity of an algebra, if given
# - treat (anti)commutativity of given function more efficiently
# - write and use a perl interface to libgap
# - use sparse matrices
# - intercept GAP errors
#

use strict;
use warnings;
use FindBin qw($Script);
use Getopt::Long;
use IPC::Open2;

# turn warnings into errors
$SIG{__WARN__} = sub { die ($_[0]); }; 


##################### parameters for GAP invocation ############################
# one may wish to use -p option in which GAP adheres to sort of non-visual 
# protocol of getting commands and displaying output (see
# http://www.gap-system.org/mailman_archives/forum/2007/001923.html ), but
# it does not seem to work properly (as of gap 4.7.5): some (simple) commands 
# apparently misbehave in that mode

use constant GAP_DIR => '/usr/local/gap';
use constant GAP_INVOCATION => 
  (GAP_DIR() . '/gap.exe', '-l', GAP_DIR(), '-b', '-q', '-n', '-E', '-T'); 

#use constant GAP_DIR => '/usr/local/gap-4.10.2';
#use constant GAP_INVOCATION => 
#  (GAP_DIR() . '/gap', '-l', GAP_DIR(), '-b', '-q', '-n', '-E', '-T');
################################################################################


use constant LEFT  => 0;
use constant RIGHT => 1;

sub _add_matrix_elem ($$$);
#sub _add_matrix_elem1 ($$$$$);
sub _add_matrix_elem2 ($$$$);
#sub _add_vector_elem ($$$);
sub _combinations ($$);
sub _debug_gap ($);
sub _debug_tmp ($);
sub _debug_trans ($);
sub _exit ($);
sub _eval ($);
sub _mult ($$);
sub _println ($);
sub _read_gap ();
sub _read_gap_notimeout ();
sub _str2ind ($);
sub _write_gap ($);
sub _writeln_gap ($);
sub _write_silent_gap ($);
sub _writeln_silent_gap ($);
sub _N (@);


##################### COMMAND LINE #############################################
my ($eqfile, $modfile, $gapfile, $logfile);
my @debug = ();
GetOptions (
			'eq=s'     => \$eqfile,
			#'modulo=s' => \$modfile,   
			'alg=s'    => \$gapfile,
			'log=s'    => \$logfile,
			'debug=s'  => \@debug,
			'version'  => sub
			{
			  print (STDERR "$Script version 0.4.3\n");
			  exit;
			},
			'help'  => sub
			{

#Usage:   $Script --eq=<file> [--modulo=<file>] --alg=<file> [--debug=gap|trans]
#--modulo=<file> <file> contains equation(s) modulo which it should be solved 
#                (e.g., like cocycles modulo coboundaries)
#          $Script --eq=2cocycle --modulo=2cobound --alg=W1.gap          

			  print <<EOF
Usage:   $Script --eq=<file> --alg=<file> [--debug=gap|trans] [--log=<file>]
         $Script --help|--version
Example: $Script --eq=der --alg=W1.gap --log=gap.log --debug=trans

--eq=<file>     <file> contains equation(s) to be solved;
--alg=<file>    <file> contains GAP code with definition of an algebra and a 
                module over it;
--debug         print debug messages: 
                gap   - related to interaction with GAP
                trans - performing perl transformations of algebraic 
                expressions;
                this option can be repeated;
--log=<file>    log GAP session to <file>;
--help          print this help;
--version       print version.

For further info, see https://www1.osu.cz/~zusmanovich/soft/lineq/
EOF
				;
			  exit;
			}
		   );

_exit ('--eq option is mandatory')  if (! defined ($eqfile));
_exit ('--alg option is mandatory') if (! defined ($gapfile));

my $DEBUG_GAP   = 0;
my $DEBUG_TRANS = 0;
for (@debug)
  {
	if ($_ eq 'gap')
	  {
		$DEBUG_GAP = 1;
	  }
	elsif ($_ eq 'trans')
	  {
		$DEBUG_TRANS = 1;
	  }
	else
	  {
		_exit ('--debug option assumes one of two values: \'gap\' or \'trans\'');
	  }
  }
################################################################################



##################### START GAP ################################################
# start GAP and perform initial calculations

my $pid_gap = open2 (*GAP_OUTPUT, *GAP_INPUT, GAP_INVOCATION);
die ('error starting GAP') if (! defined ($pid_gap));
_debug_gap ('GAP started');

if (defined ($logfile))
  {
	_writeln_gap ("LogTo (\"${logfile}\")");
	_debug_gap ("logging to $logfile");
  }

# feed the content of the second argument to GAP
open (FILE, "< $gapfile") || die ("error opening $gapfile");
my $last_line;
for (<FILE>)
  {
	s/#.*$//g;        # ignore comments
	next if /^\s*$/;  # skip blank lines
	$last_line = $_;
	#_write_gap ($last_line);
	_write_silent_gap ($last_line);
  }
close (FILE);

# the last (sans blanks and comments) line at this file should have a form 
# M := ...
# where M is a module in which the unknown function takes the values
my $module_var;
if ($last_line =~ /^\s*([A-Za-z][A-Za-z0-9]*)\s*:=/)
  {
	$module_var = $1;
  }
else
  {
	_exit ('cannot determine GAP variable for module');
  }

## the number of times we should read the GAP input
## could be changed based on the input
#_read_gap();

_writeln_silent_gap ("A := LeftActingAlgebra (${module_var})");
_writeln_silent_gap ('b := CanonicalBasis (A)');
_writeln_gap ('Dimension (A)');

# here we reading from GAP for the first time; do not set timeout, as if not
# in the cached memory, the start of GAP can take a considerable amount of time 
my $dim_alg = _read_gap_notimeout();

# we do not use CanonicalBasis() here, as some modules do not have 
# CanonicalBasis
_writeln_silent_gap ("m := Basis (${module_var})");
_writeln_gap ("Dimension (${module_var})");
my $dim_mod = _read_gap();  # dimension of the module
################################################################################



##################### EQUATIONS: VARIABLES, ARITY, ETC. ########################
# read a file containing equations to be solved
# one equation should occupy one line; blank lines and comments are ignored
open (FILE, "< $eqfile") || die ("error opening $eqfile");
my @lhs = ();  # list of left-hand sides of equations to be solved
_debug_trans ("solving equation(s):\n");
for (<FILE>)
  {
	s/#.*$//g;        # ignore comments
	s/\s*//g;         # trim white space
	next if /^\s*$/;  # skip blank lines
	push (@lhs, $_);
	_debug_trans ("$_ = 0\n");
  }
close (FILE); 

# values for $status
use constant NOT_INSIDE_FUNCTION => 0;
use constant AFTER_FUNCTION      => 1;
use constant INSIDE_FUNCTION     => 2;
use constant INSIDE_FUNCTION_READY_TO_CLOSE => 3;

my %LHS = ();  # left-hand side of the equation => number of variables
my $N = 0;
my ($lhs, $lhs_new, $map_symb, $arity, $a, @a, %vars, $num_vars, $status);
for $lhs (@lhs)
  {
	# determine which variables (elements of an algebra) are used (x,y,z, ...)
	#($a = $lhs) =~ s/[^a-z]//g;
	($a = $lhs) =~ s/[^ac-z]//g;
	%vars = ();
	for (split (//, $a))
	  {
		$vars{$_} = 1;
	  }

	# determine the arity of the map
	# pasha sep 3 2019: non-greedy matching, only letters, comma and parentheses
	$lhs =~ /([A-Z])\(([a-z,\*\(\)]+)\)/;
	if (defined ($map_symb))
	  {
		# be sure the same symbol is used for the map in all equations
		_exit 
		  ("ERROR: different symbols $map_symb and $1 are used in different equations")
			if ($map_symb ne $1);
	  }
	else
	  {
		$map_symb = $1;
	  }
	@a = split (/,/, $2);
	if (defined ($arity))
	  {
		_exit ('ERROR: different arities in different equations') 
		  if ($arity ne scalar (@a));
	  }
	else
	  {
		$arity = scalar (@a);
	  }

	# replace parentheses in F(...) by curly braces: F{...} 
	# (to facilitate subsequent parsing)
	# pasha sep 6 2019: switched to character-by-characer processing instead
	# of regex replacement
	#$lhs =~ s/${map_symb}\(([^${map_symb}]+?)\)/${map_symb}{$1}/g;
	#$lhs =~ s/${map_symb}\(([^${map_symb}]+)\)/${map_symb}{$1}/g;
	
	$lhs_new = '';
	my $comma_count = 0;
	$status = NOT_INSIDE_FUNCTION;
	for (split (//, $lhs))
	  {
		if ($status == NOT_INSIDE_FUNCTION)
		  {
			if ($_ eq $map_symb)
			  {
				$status = AFTER_FUNCTION;
				$comma_count = 0;
			  }
			$lhs_new .= $_;
			next;
		  }
		if ($status == AFTER_FUNCTION)
		  {
			if ($_ eq '(')
			  {
				$status = INSIDE_FUNCTION;
				$lhs_new .= '{';
				next;
			  }
			else
			  {
				_exit ("malformed equation $lhs");
			  }
		  }
		if ($status == INSIDE_FUNCTION)
		  {
			if ($_ eq ',')
			  {
				$comma_count++;
			  }
			if ($comma_count == $arity - 1)
			  {
				$status = INSIDE_FUNCTION_READY_TO_CLOSE;
			  }
			$lhs_new .= $_;
			next;
		  }
		if ($status == INSIDE_FUNCTION_READY_TO_CLOSE)
		  {
			if ($_ eq ',')
			  {
				_exit ("wrong number of commas in $lhs"); 
			  }
			if ($_ eq ')')
			  {
				$lhs_new .= '}';
				$status = NOT_INSIDE_FUNCTION;
			  }
			else
			  {
				$lhs_new .= $_;
			  }
			next;
		  }
	  } # end of loop thru symbols of $lhs

	$lhs = $lhs_new;

	# replace each variable in the left-hand side by an element of the basis 
	# of an algebra
	# (x -> e_i0, y -> e_i1, z -> e_i2, ...)
	$num_vars = 0;
	for (keys %vars)
	  {
		$lhs =~ s/$_/e_i${num_vars}/g;
		$num_vars++;
	  }
	$N += $dim_alg ** $num_vars;  # used later to compute total number of eqs
	$LHS{$lhs} = $num_vars;
  } # end of loop thru equations
################################################################################



##################### REGEXPS ##################################################
# some regexps used below
# TODO: compile them beforehand?
my $re_int = '[1-9][0-9]*';                    # a positive integer
my $re_rat = "(?:$re_int|${re_int}/$re_int)";  # a positive rational
my $re_b   = "b\[${re_int}\]";                 # b[i]

# determine the base field
_writeln_gap ("K := LeftActingDomain (A)");
my $base_field = _read_gap();

my $re_coeff;   # element of the base field, surrounded by brackets
if ($base_field =~ /^GF\(/)
  {
	# finite field
	$re_coeff = "\\((?:(?:|\-)${re_rat}|Z\\(${re_int}\\)(?:|\\^(?:${re_int}|0)))\\)";
  }
else
  {
	# rationals
	$re_coeff = "\\((?:|\-)${re_rat}\\)";
  }
################################################################################



##################### INITIALIZE MATRIX ########################################
# Construct matrix for the homogeneous linear system.
# There are ($dim_alg^$arity) * $dim_mod unknowns and 
# ($dim_alg^$num_vars_1 + ... + $dim_alg^$num_vars_n) * $dim_mod equations,
# where $num_vars_1, ..., $num_vars_n are number of variables in each
# of the equations.
#
# Unknowns are coefficients x_{i1,i2,...,i_arity,j} in decompositions
# F(b_i1, b_i2, ..., b_i_arity) = 
#       \sum_{j=1}^dim_mod x_{i1,i2,...,i_arity,j} m_j
#
# They are arranged as one row in the matrix in lexicographic order:
# x_{11...11}, x_{11...12}, ..., x_{11...1,dim_mod}, 
# x_{11...21}, x_{11...22}, ..., x_{11...2,dim_mod}, 
# ..., 
# x_{dim_alg,dim_alg,...,dim_alg,dim_mod},
# so the element x_{i1,i2,...,i_arity,j} is on the place with the number
# N(i1,...,i_arity)*dim_mod + j, where
#
# N(i1,...,i_arity)
# = (i1 - 1)          * dim_alg^(arity-1)
# + (i2 - 1)          * dim_alg^(arity-2)
# + ... 
# + (i_{arity-1} - 1) * dim_alg
# + (i_arity - 1) 
#
# Equations coming from substituting into each $lhs every
# possible combination of elements of the basis of the algebra, and
# equating coefficients by the basis of the module to zero.

## In the case the module is the algebra itself, we may also have nonhomogenous 
## components, i.e. elements of the basis of the algebra.
## If an element coeff*b[i] is encountered in equation number j, then the vector
## Vec in the nonhomogeneous linear system X Mat = Vec acquires -coeff at ? place

my $num_unknowns = ($dim_alg ** ${arity}) * $dim_mod;
my $num_eqs      = $N * $dim_mod;
_writeln_silent_gap ("Mat := NullMat (${num_unknowns}, ${num_eqs}, K)");

## NONHOMOG
#_writeln_silent_gap ("Vec := ListWithIdenticalEntries (${num_eqs}, Zero(K))");

################################################################################



##################### LOOP OVER EQUATIONS ######################################
my $eq = 0;  # total counter for equations
my ($tmp, $summand, $c);
for $lhs (keys %LHS)
  {
	_debug_trans ("processing equation ${lhs}\n");

	# loop over all combinations of basic elements occuring in $lhs
	for $c (_combinations ($LHS{$lhs}, $dim_alg))
	  {
		$eq++;
		$tmp = $lhs;
		for ($tmp)
		  {
			# substitute in current $lhs the corresponding basis elements
			s/e_i([0-9]+)/b[$c->[$1]]/g;
			_debug_tmp ('SUBS');

			# the things below seem to be a suitable job for simplification 
			# facilities of some general-purpose computer algebra system, but 
			# neither Maxima nor Reduce seem to be up to the task - 
			# the approprite simplify routines for multilinear functions handle,
			# strangely, the first argument only and with weird side-effects;
			# Maple, on the other hand, is seemingly capable of doing that.

			# replace all occurencies of products of b[i]'s by their product in 
			# the algebra
			# pasha sep 7 2019: allow multiple products with parenthesis --
			# like (b[1]*b[2])*b[3]

			# the next line is from version 0.4.2, introduced on aug 18 2020, 
			# and it breaks some things (e.g., hom-lie)
			#s/([b\[\]0-9\(\)]+\*[b\[\]0-9\(\)\*]+)/_eval ($1)/eg;
			s/([b\[\]0-9\(\)]+\*[b\[\]0-9\(\)]+)/_eval ($1)/eg;

			_debug_tmp ('PROD');

			while (
				   # replace product with zero by zero, i.e.:
				   # (0)*b[n] = 0, b[n]*0 = 0
				   (s/\(0\)\*${re_b}/0/g) ||
				   (s/${re_b}\*0/0/g)     ||

				   # replace functions with one of arguments is zero, by zero 
				   # (i.e. F(0,*) = 0, etc.)
				   (s/${map_symb}\{(?:|[^\{\}]+?,)0(?:|,[^\{\}]+?)\}/0/g) ||

				   # replace terms which are multiplies of zero, by zero;
				   # i.e.: (0)*F(...) = 0, 0*F(...), and F(...)*(0)
				   #(s/\(0\)\*${map_symb}\{[^\{\}]+?\}/0/g) ||
				   (s/(\(0\)|0)\*${map_symb}\{[^\{\}]+?\}/0/g) ||
				   (s/${map_symb}\{[^\{\}]+?\}\*\(0\)/0/g) ||

				   # get rid of summand which are zeros
				   (s/(^0|\+0|-0)//g)
				  )
			  {
				_debug_tmp ('ZERO');
			  }

			# (coef*b[i] + coef*b[j] + ...)*F = coef*b[i]*F + coef*b[j]*F + ...
			while (
				   s/\(([\(\)\-\*\^bZ\[\]0-9]+?)\+([\(\)\-\+\*\^bZ\[\]0-9]+?)\)\*(${map_symb}\{[^\{\}]+?\})/($1)*$3+($2)*$3/g
				  )
			  {
				_debug_tmp ('DIST');
			  }

			# remove redundant parentheses:
			# (b[i]) -> b[i]
			if (s/\((${re_b})\)/$1/g)
			  {
				_debug_tmp ('PARN');
			  }
			#s/\(($re_coeff)\)/$1/g;
			if (s/\((${re_coeff}\*${re_b})\)/$1/g)
			  {
				_debug_tmp ('PARN');
			  }

			# replace coeff* by (coeff)* and
			# -coeff* by +(-coeff)*
			s/\-(${re_rat})\*/+(-$1)*/g;
			s/(${re_rat})\*/($1)*/g;
			s/-\((${re_rat})\)/+(-$1)/g;

			# replace '-' by (-1)*
			if (s/\-(${map_symb}|b)/+(-1)*$1/g)
			  {
				_debug_tmp ('MNUS');
			  }
			
			# expand by linearity
			# F(...,a+b) = F(...,a) + F(...,b)
			# non-greedy matching
			while (
				   #s/${map_symb}\{(|.*?\,)(|${re_coeff}\*)(${re_b})\+([bZ\[\]\(\)\*\-\+\^0-9]+?)\}/${map_symb}{$1$2$3}+${map_symb}{$1$4}/g
				   s/(|${re_coeff}\*)${map_symb}\{(|.*?\,)(|${re_coeff}\*)(${re_b})\+([bZ\[\]\(\)\*\-\+\^0-9]+?)\}/$1${map_symb}{$2$3$4}+$1${map_symb}{$2$5}/g
				  )
			  {
				#print (STDERR "YOYO1 \$1=$1 \$2=$2 \$3=$3 \$4=$4 \n");

				#D{b[13]+b[14]}
				#$1= $2= $3=b[13] $4=b[14] $5=} 
				
				_debug_tmp ('LIN0');
				#exit;
			  }
			
			# a variant of the previous
			# F(...,a+b,...) = F(...,a,...) + F(...,b,...)
			while (
				   s/(|${re_coeff}\*)${map_symb}\{(|.*?\,)(|${re_coeff}\*)(${re_b})\+(.*?)\,(.*?)\}/$1${map_symb}{$2$3$4,$6}+$1${map_symb}{$2$5,$6}/g
				  )
			  {
				#print (STDERR "YOYO2 \$1=$1 \$2=$2 \$3=$3 \$4=$4 \$5=$5 \n");

				#F{(-1)*b[13]+(-1)*b[14],b[1]}
				#$1= $2=(-1)* $3=b[13] $4=(-1)*b[14] $5=,
				
				_debug_tmp ('LIN1');
			  }

			# F(...,la,...) = lF(...,a,...)
			while 
			  (
			   s/${map_symb}\{(|[^\{\}]+?,)(${re_coeff})\*([^\{\},]+?)(|,[^\{\}]+?)\}/$2*${map_symb}{$1$3$4}/g
			  )
			  {
				_debug_tmp ('LIN2');
			  }

			# evaluate products of coefficents
			if (s/(${re_coeff})\*(${re_coeff})/'(' . _mult ($1, $2) . ')'/eg)
			  {
				_debug_tmp ('COEF');
			  }
		  } # end of transformations of $tmp
		
		_debug_trans ("\n");

		# if $tmp evaluates entirely to zero, at this point it will be empty
		next if ($tmp =~ /^\s*$/);
		
		# at this point $tmp should consist of the sum of expressions of the 
	    # form:
		#
		# coeff * b[k1] * b[k2] * ... * F(b[i],b[j],...) * b[s1] * b[s2] * ...
		#
		# with some bracketing; coeff may be absent, and external b[k]'s may be
		# absent -- from the left, from the right, or from both sides, i.e., 
		# having just
		#
		# coeff * F(b[i],b[j],...) or even just F(b[i],b[j],...)
		
		# split the sum into summands
		@a = split (/\+/, $tmp);
		
		# if the first summand is with the plus sign, 
		# the first matching pattern will be empty
		shift (@a) if ($a[0] eq '');
		
		for $summand (@a)  # loop thru summands in a current expression
		  {
			for ($summand)
			  {
				## NONHOMOG
				#/^(${re_coeff})\*([^\{\}]+?)$/ && do
				#  {
				#	# coeff * b[i]
				#	print "processing summand $summand : $1 $2 \n";
				#	exit;
                #
				#	# Vec[($eq - 1)*$dim_mod + j] = -coeff
				#	_add_vector_elem ($1, $eq, $2);
                #
				#	last;
				#  };

				/^${map_symb}\{([^\{\}]+?)\}$/ && do
				  {
					# F{b[i],...}
					#print "processing summand $summand : $1 \n";
					
					# same as the next one, with coeff = 1
					_add_matrix_elem ($1, $eq, 1);
					last;
				  };
				
				/^(${re_coeff})\*${map_symb}\{([^\{\}]+?)\}$/ && do
				  {
					# coeff * F{b[i],...}
					#print "processing summand $summand : $1 $2 \n";
					
					# Let F(b_i1, b_i2, ..., b_i_arity) = 
					#       \sum_{j=1}^dim_mod x_{i1,i2,...,i_arity,j} m_j
					# where {m_j} - basis of the module M.
					# Set
					# Mat[N(i1,...,i_arity)*dim_mod) + j][($eq - 1)*$dim_mod + j]
					#        += coeff
					# where N(i1,...,i_arity) is the number as above.
					
					_add_matrix_elem ($2, $eq, $1);
					last;
				  };
				

				# pasha sep 22-23,dec 8-9 2019: 
				# incorporated possible scalar coefficient; 
				# added ^ (module action)
				#/^(${re_coeff}\*|)([b\[\]0-9\*\(\)]*)${map_symb}\{([^\{\}]+?)\}([b\[\]0-9\*\(\)]*)$/ && do
				($base_field =~ /^GF\(/ ? 
				/^([b\[\]0-9\*\^\(\)\-Z]*)${map_symb}\{([^\{\}]+?)\}([b\[\]0-9\*\^\(\)]*)$/ 
				:
				/^([b\[\]0-9\*\^\(\)\-\/]*)${map_symb}\{([^\{\}]+?)\}([b\[\]0-9\*\^\(\)]*)$/ 
				) && do
				  {
					# expressions of the form
					#
					# coeff * b[k1] * b[k2] * ... * F(...) * b[s1] * b[s2] * ...
					#
					# where instead of * might stand ^, whether appropriate,
					# with appropriate bracketing -- i.e., multiplication/action
					# on F(...) from the left and from the right
					
					# Let F(b_i1, b_i2, ..., b_i_arity) = 
					#       \sum_{j=1}^dim_mod x_{i1,i2,...,i_arity,j} m_j
					# and 
					# coeff * b[k1] * b[k2] * ... * F(...) * b[s1] * b[s2] * ...
                    #        = \sum_{s=1}^{dim_mod} l_ks m_s
					# where at the left-hand side is the experssion above, and
					# {m_j} - basis of the module M.
					# Then 
					# coeff * b[k1] * b[k2] * ... * F(...) * b[s1] * b[s2] * ...
					# = \sum_{s=1}^{dim_mod} (\sum{j=1}^dim_mod coeff*l_ks)
					#        x_{i1,i2,...,i_arity,j} m_s

					#print $1, "\n";
					#print $2, "\n";
					#print $3, "\n";
					#exit;

					_add_matrix_elem2 ($1, $2, $3, $eq);
					last;
				  };

				_exit ("ERROR: wrong summand $summand in expression $tmp");
			  }
		  } # end of loop over summands
	  } # end of loop over all combinations of basic elements
  } # end of loop over equations
################################################################################



##################### SOLVE THE LINEAR SYSTEM ##################################
# TODO: use other system/package for solving linear systems, more efficient 
# than GAP?

#_writeln_silent_gap ("
#PrintTo (\"${outfile}\", NullspaceMatDestructive (Mat), \"\\n\")
#");

print (STDERR "finished compose linear system, solving it ...\n");
_writeln_gap ('Print (NullspaceMatDestructive (Mat), "\n")');

## NONHOMOG
##_writeln_gap ('Print (SolutionMatDestructive (Mat, Vec), "\n")');
##_read_gap();


################################################################################



##################### PRINT THE ANSWER #########################################
# get the GAP output;
# as GAP insists on formatting output and breaking it line-by-line,
# and the standard things like checking whether stream is ready for read/write
# do not work with it, we resort here on an ugly, output-specific, way to 
# determine when GAP output is finished
my $str = '';
my $line;
my $first_time = 1;
while (1)
  {
	# for the first time, do not set timeout for reading GAP, as it solves
	# linear system which may takes minutes and more!
	if ($first_time)
	  {
		$line = _read_gap_notimeout();
		$first_time = 0;
	  }
	else
	  {
		$line = _read_gap();
	  }
	$line =~ s/\s*//g;  # remove white space
	if ($line =~ /\[\]/)
	  {
		# the solution of linear system is zero (should be on the first line)
		if ($str ne '')
		  {
			_exit ('error in GAP output');
		  }
		print ("dimension: 0\n");
		_exit (undef);
	  }
  	$str .= $line;
	last if ($line =~ /\]\]/)
  }

# write each row in the solution matrix as the map on the algebra, 
# in terms of the basis
for ($str)
  {
	s/^\[\[//;
	s/\]\]$//;
  }
@a = split (/\],\[/, $str);
print ('dimension: ' . scalar (@a) . "\n\n");
my @basis_combinations = _combinations ($arity, $dim_alg);

my ($j, $index, $val, @b);
my $count = 0;
for (@a)
  {
	$count++;
	print ("#${count}:\n");

	@b = split (/,/);
	for $c (@basis_combinations)
	  {
		print (join (' x ', map {'e['. $_ . ']'} @{$c}) . ' -> ');
		$str = '';
		for $j (1 .. $dim_mod)
		  {
			$index = _N (@{$c}) * $dim_mod + $j;
			$val = $b[$index - 1];
			$str .= "${val}*m[${j}] + " if (($val ne '0') && ($val !~ /^0\*Z/));
		  }
		$str =~ s/\s*\+\s*$//;
		print (($str eq '' ? 0 : $str), "\n");
	  }
	print ("\n");
  } # end of loop over solutions of the linear system
################################################################################


_exit (undef);



##################### FUNCTIONS ################################################
# evaluate expression which supposed to be operation on basic elements --
# e.g. (b[1]*b[2])*b[3] -- in GAP, return the result in the form of linear
# combinations of b[n]'s
sub _eval ($)
  {
	# check that parentheses are balanced
	my $count_left  = $_[0] =~ tr/(//;
	my $count_right = $_[0] =~ tr/)//;
	return ($_[0]) if ($count_left != $count_right);

	_writeln_gap ("Coefficients (b, $_[0])");
	my $result = '';
	my $line;
	while (1)
	  {
		$line = _read_gap();
		$result .= $line;
		last if ($line =~ /\w*\]/);
	  }

	# $result is a GAP array (comma-delimited, surrounded by square brackets)
    # containing coefficients of the algebra element in the basis b
	for ($result)
	  {
		s/^\[\s*//;
		s/\s*\]$//;
	  }
	my $i = 0;
	my @a = ();
	for (split (/\s*,\s*/, $result))
	  {
		$i++;
		next if (/^0/);  # zero coefficient
		push (@a, "(${_})*b[${i}]");
	  }

	_exit 
	  ("ERROR: the number of coefficient in the result of evaluation of $_[0] does not coincide with the dimension of algebra")
	  if ($i != $dim_alg);
	return ((scalar (@a) > 0) ? join ('+', @a) : '0');
  }


# multiply in GAP two elements of the ground field
sub _mult ($$)
  {
	_writeln_gap ($_[0] . '*' . $_[1]);
	return (_read_gap());
  }


# add elements to GAP matrix of the homegenous system, case coeff*F(b[i],...)
# args:
# 0 - a string b[i],... 
#     (a number, depending on the arity, of b[i]'s separated by commas)
# 1 - current sequential number of basis elements combination ("equation 
#     number")
# 2 - coeff
sub _add_matrix_elem ($$$)
  {
	my $I = _N (_str2ind ($_[0])) * $dim_mod;
	my $J = ($_[1] - 1) * $dim_mod;
	my $elem_to_add = $_[2];
	_writeln_silent_gap ("
for j in [1 .. ${dim_mod}] do 
    Mat[${I}+j][${J}+j] := Mat[${I}+j][${J}+j] + ${elem_to_add};
od
");
  }


# pasha dec 9 2019: removed, no longer used
## ditto, case 
## coeff * b[k] * F(b[i],...) (left action) or
## coeff * F(b[i],...) * b[k] (right action)
## args: 
## 0 - k (index in b[k])
## 1 - a string b[i],... 
##     (a number, depending on the arity, of b[i]'s separated by commas)
## 2 - current sequential number of basis elements combination ("equation 
##     number")
## 3 - coeff
## 4 - whether action of b[k] is from the left  (LEFT),  i.e. b[k]*F(...), 
##                            or from the right (RIGHT), i.e. F(...)*b[k]
#sub _add_matrix_elem1 ($$$$$)
#  {
#	my $k = $_[0];
#	my $I = _N (_str2ind ($_[1])) * $dim_mod;
#	my $J = ($_[2] - 1) * $dim_mod;
#	my $coeff = $_[3];
#	if ($_[4] == LEFT)
#	  {
#		_writeln_silent_gap ("
#for j in [1 .. ${dim_mod}] do
#   	lambda := Coefficients (m, b[${k}]^m[j]);
#	for s in [1 .. ${dim_mod}] do
#		Mat[${I}+j][${J}+s] := Mat[${I}+j][${J}+s] + ${coeff}*lambda[s];
#	od;
#od
#");
#	  }
#	elsif ($_[4] == RIGHT)
#	  {
#		_writeln_silent_gap ("
#for j in [1 .. ${dim_mod}] do
#   	lambda := Coefficients (m, m[j]^b[${k}]);
#	for s in [1 .. ${dim_mod}] do
#		Mat[${I}+j][${J}+s] := Mat[${I}+j][${J}+s] + ${coeff}*lambda[s];
#	od;
#od
#");
#	  }
#	else
#	  {
#		_exit ('wrong parameter ' . $_[4]);
#	  }
#  } # end of _add_matrix_elem1()




# add elements to GAP matrix of the homegenous system in the most general case:
# b[k1] * b[k2] * ... * F(...) * b[s1] * b[s2] * ...
# with appropriate bracketing
# args:
# 0 - the left part b[k1] * b[k2] * ... *, with brackets
# 1 - a string b[i],... 
#     (a number, depending on the arity, of b[i]'s separated by commas)
# 2 - the right part * b[s1] * b[s2] * ..., with brackets
# 3 - current sequential number of basis elements combination ("equation 
#     number")
sub _add_matrix_elem2 ($$$$)
  {
	my $I = _N (_str2ind ($_[1])) * $dim_mod;
	my $J = ($_[3] - 1) * $dim_mod;
	my $coeff = 1;

	# form a GAP string of b[k]'s acting on m[j] from left and right, as in
	# $_[0] and $_[2]
	my $action = $_[0] . 'm[j]' . $_[2];

	_writeln_silent_gap ("
for j in [1 .. ${dim_mod}] do
   	lambda := Coefficients (m, ${action});
	for s in [1 .. ${dim_mod}] do
		Mat[${I}+j][${J}+s] := Mat[${I}+j][${J}+s] + ${coeff}*lambda[s];
	od;
od
");
  } # end of _add_matrix_elem2()




## NONHOMOGEN
## add element to the vector in the linear system (nonhomogeneous part)
## args: 
## 0 - a string of the form "b[i]"  
## 1 - current sequential number of basis elements combination ("equation 
##     number")
## 2 - coeff
#sub _add_vector_elem ($$$)
#  {
#	my $J = ($_[1] - 1) * $dim_mod;
#	#Vec[$J] := -coeff;
#
#  } # end of _add_vector_elem()




sub _N (@)
  {
	my $N = 0;
	for (@_)
	  {
		$N = $N*$dim_alg + $_ - 1;
	  }
	return ($N);
  }


# converts a string "b[i],..." to array of indices (i, ...)
sub _str2ind ($)
  {	
	my @arr = ();
	for my $b (split (/,/, $_[0]))
	  {
		# extract i from b[i]
		for ($b)
		  {
			s/^b\[//;
			s/\]$//;
		  }
		push (@arr, $b);
	  }
	_exit ("ERROR: string of arguments $_[0] is not of arity $arity")
	  if (scalar (@arr) != $arity);
	return (@arr);
  }


# produces all combinations of n (0th argument) symbols of size k (1st argument)
# (e.g., n=3 k=2 -> 11 12 13 21 22 23 31 32 33)
sub _combinations ($$)
  {
	my ($n, $size) = @_;
	my @combinations = ();
	my @a = ();
	my $loop = 'my (';
	my $i;
	for $i (0 .. $n - 1)
	  {
		$loop .= "\$i${i}, ";
	  }
	$loop .= '); ';
	for $i (0 .. $n - 1)
	  {
		$loop .= "for \$i$i (1 .. ${size}) { push (\@a, \$i$i); ";
	  }
	$loop .= 'push (@combinations, [@a]); ';
	for $i (0 .. $n - 1)
	  {
		$loop .= ' pop (@a); } ';
	  }
	eval ($loop);
	_exit ("error evaluating nested loop: ${loop}: $@") if ($@);
	return (@combinations);
  }


# write line to GAP input as is
sub _write_gap ($)
  {	
	_debug_gap ('writing to GAP: ' . $_[0]);
	print (GAP_INPUT $_[0]);
  }

# write line to GAP, supplying the end of line 
# (normally used when expecting output from it)
sub _writeln_gap ($)
  {
	_write_gap ($_[0] . ";\n");
  }

# write line to GAP, supressing output
sub _writeln_silent_gap ($)
  {
	_write_gap ($_[0] . ";;\n");
  }

# write line to GAP, supressing output only if this is end of a command
# (i.e. ends with ';'); otherwise, write as is
sub _write_silent_gap ($)
  {
	$_[0] =~ s/\s*$//;
	$_[0] .= ";\n" if ($_[0] =~ /;$/);
 	_write_gap ($_[0]);
  }

sub _read_gap ()
  {
	my $line;
	# allow no more than 5 second for reading from GAP; 
	# from Perl Cookbook, 16.21
	eval 
	  { 
		local $SIG{ALRM} = sub { die ('__ALARM_CLOCK_RESTART') };
		alarm 5;              # schedule alarm in 5 seconds 
		$line = <GAP_OUTPUT>;
		chomp ($line);        # remove trailing newline
		$line =~ s/^.*\r//;   # remove some cruft GAP put sometimes in front
		_debug_gap ("reading from GAP: $line");
		alarm 0;
	  };
	alarm 0;                  # race condition protection
	_exit ('ERROR: timeout when reading GAP') 
	  if ($@ && ($@ =~ /__ALARM_CLOCK_RESTART/));
	return ($line);
  }


sub _read_gap_notimeout ()
  {
	my $line;
	$line = <GAP_OUTPUT>;
	chomp ($line);     
	$line =~ s/^.*\r//;
	_debug_gap ("reading from GAP: $line");
	return ($line);
  }


sub _exit ($)
  {
	# stop GAP
	if (defined ($pid_gap))
	  {
		close (GAP_INPUT);
		close (GAP_OUTPUT);
		waitpid ($pid_gap, 0);
		_debug_gap ('GAP stopped');
	  }
	_println ($_[0]) if (defined ($_[0]));
	exit;
  }

sub _println ($)
  {
	print (STDERR $_[0] . "\n");
  }

sub _debug_gap ($)
  {
	_println ($_[0]) if ($DEBUG_GAP == 1);
  }

sub _debug_trans ($)
  {
	print (STDERR $_[0]) if ($DEBUG_TRANS == 1);
  }

sub _debug_tmp ($)
  {
	_debug_trans ('(' . $_[0] . ") $tmp \n") if ($tmp !~ /^\s*$/);
  }
################################################################################


__END__
