#
# file:    poisson.gap
#
# purpose: Constructs Poisson algebras and their deformations. 
#          This is accompanying code to the paper
#          "Melikyan algebra is a deformation of a Poisson algebra".
#
#          Routines defined here:
#
#          DividedPowersAlgebra
#          MapOtimesOne
#          MultipleOfAlgebraStructure 
#          OneOtimesMap
#          PoissonAlgebra
#          RightMult
#          Shift
#          SumOfAlgebraStructures
#          SumOfTwoAlgebraStructures
#          TensorProductOfAlgebras
#          TensorProductOfTwoAlgebras
#
# created: Pasha Zusmanovich <pasha.zusmanovich@gmail.com> Jan 7 2014
# 
# revision history:
# Sep 25 2014,Oct 25 2015,Oct 20 2021  comments
#


# returns a divided powers algebra O_1(n) over a field K:
# this is a commutative associative algebra with the basis {x^i | 0 <= i < p^n}
# and multiplication x^i x^j = \binom{i+j}{j} x^{i+j}
DividedPowersAlgebra := function (K, n)
	local p, pn, i, j, T, O;

	p := Characteristic(K);
	if (p = 0) then
		Error ("the ground field supposed to be of positive characteristic");
	fi;
	pn := p^n;
	T := EmptySCTable (pn, Zero(K), "symmetric");

	# the basis element x^i (0 <= i <= p^n-1) goes to the place (i+1)
	for i in [0 .. pn-1] do
		for j in [0 .. i] do
			if (i+j < pn) then
				# this corresponds to x^i x^j = \binom{i+j}{j} x^{i+j}
				SetEntrySCTable (T, i+1, j+1, [Binomial (i+j,j), i+j+1]);
			fi;
		od;
	od;

	O := AlgebraByStructureConstants (K, T);
	SetIsAssociative (O, true);
	return (O);
end;




# returns an algebra structure which is a multiple of a given one
MultipleOfAlgebraStructure := function (A, l)
	local K, n, i, j, limit_for_j, k, T, sc, list, B;	

	K := LeftActingDomain(A);
	n := Dimension (A);

	# structure constants of the resulting algebra
	if (IsAnticommutative (A)) then
		T := EmptySCTable (n, Zero(K), "antisymmetric");
	elif (IsCommutative (A)) then	
		T := EmptySCTable (n, Zero(K), "symmetric");
	else
		T := EmptySCTable (n, Zero(K));
	fi;

	sc := StructureConstantsTable (Basis (A));

	limit_for_j := n;
	for i in [1 .. n] do
		if (IsCommutative (A)) then
			limit_for_j := i;
		elif (IsAnticommutative (A)) then
			limit_for_j := i-1;
		fi;
		for j in [1 .. limit_for_j] do
			list := [];
			for k in [1 .. Length (sc[i][j][1])] do
				Add (list, sc[i][j][2][k] * l);
				Add (list, sc[i][j][1][k]);
			od;
			SetEntrySCTable (T, i, j, list);
		od;
	od;

	B := AlgebraByStructureConstants (K, T);
	if (IsAssociative (A)) then
		SetIsAssociative (B, true);
	elif (IsLieAlgebra (A)) then
		SetIsLieAlgebra (B, true);
	fi;

	return (B);
end;




# for an algebra A, and two maps D,F: A->A, builds an algebra on the vector space
# A with multiplication a*b = D(a)F(b) - F(a)D(b)
PoissonAlgebra := function (A, D, F)
	local K, B, T, i, j, limit_for_j, list, cf, n;

	# IsLinearMapping() is not implemented as of GAP 4.7.2
	if ((not IsLeftModuleHomomorphism (D)) or (not IsLeftModuleHomomorphism (F)))
	then
		Error ("both arguments supposed to be linear maps");
	fi;
	if ((Source (D) <> A) or (Source (F) <> A)) then
		Error ("sources of both arguments should coincide with an algebra");
	fi;

	K := LeftActingDomain (A);
	B := Basis (A);

	# structure constants of the resulting algebra
	if (IsCommutative (A)) then
		T := EmptySCTable (Dimension(A), Zero(K), "antisymmetric");
	else 
		T := EmptySCTable (Dimension(A), Zero(K));
	fi;

	limit_for_j := Dimension (A);
	for i in [1 .. Dimension(A)] do
		if (IsCommutative (A)) then
			limit_for_j := i-1;
		fi;
		for j in [1 .. limit_for_j] do
			list := [];
			n := 0;
			for cf in Coefficients (B, Image (D, B[i]) * Image (F, B[j]) - 
									   Image (F, B[i]) * Image (D, B[j]))
			do
				n := n+1;
				if (cf <> Zero(K)) then
					Add (list, cf);
					Add (list, n);
				fi;
			od;
			SetEntrySCTable (T, i, j, list);
		od;
	od;

	# TODO: if D and F are commuting derivations, then set the resulting algebra
	# as a Lie algebra (probably do that if an optional argument is supplied?)

	return (AlgebraByStructureConstants (K, T));
end;




# for an algebra A and its element a, builds a linear map on A which is
# a right multiplication by a
RightMult := function (A, a)
	local B;
	B := Basis (A);
	return (LeftModuleHomomorphismByImages 
		(A, A, B, List ([1 .. Dimension(A)], i -> B[i]*a)));
end;




# for two algebra structures (A, *_A) and (B,*_B) on the same vector space, 
# return an algebra structure * = *_A + *_B
# the code resembles AntiCommAlg() (which, in its turn, is taken from 
# LieAlgebraByDomain()
SumOfTwoAlgebraStructures := function (A, B)
	local K, T, scA, scB, nullvec, aij, bij, cfs, list, i, j, k, m, n;

	K := LeftActingDomain(A);
	if (K <> LeftActingDomain(B)) then
		Error ("algebras should be defined over the same field");
	fi;
	n := Dimension(A);
	if (Dimension (B) <> n) then
		Error ("dimension of algebras should coincide");
	fi;

	# structure constants of the resulting algebra
	if (IsAnticommutative (A) and IsAnticommutative (B)) then
		T := EmptySCTable (n, Zero(K), "antisymmetric");
	elif (IsCommutative (A) and IsCommutative (B)) then	
		T := EmptySCTable (n, Zero(K), "symmetric");
	else
		T := EmptySCTable (n, Zero(K));
	fi;

	scA := StructureConstantsTable (CanonicalBasis (A));
	scB := StructureConstantsTable (CanonicalBasis (B));
	nullvec := List ([1..n], x -> Zero(K)); 

	for i in [1 .. n] do
		for j in [1 .. n] do
			cfs := ShallowCopy (nullvec);
			aij := scA[i][j];
			bij := scB[i][j];
			for m in [1 .. Length (aij[1])] do
				k := aij[1][m];
				cfs[k] := cfs[k] + aij[2][m];
			od; 
			for m in [1 .. Length (bij[1])] do
				k := bij[1][m];
				cfs[k] := cfs[k] + bij[2][m];
			od; 
			list := [];
			for m in [1..n] do
				if (cfs[m] <> Zero(K)) then
					Add (list, cfs[m]);
					Add (list, m);
				fi;
			od; 
			SetEntrySCTable (T, i, j, list);
		od;
	od;

	return (AlgebraByStructureConstants (K, T)); 
end;




SumOfAlgebraStructures := function (arg)
	local T, i;
	if (Length(arg) = 0) then
		Error ("no arguments supplied");
	fi;
	T := arg[1];
	for i in [2 .. Length(arg)] do
		T := SumOfTwoAlgebraStructures (T, arg[i]);
	od;
	return (T);
end;



# for an algebra A with a canonical basis {e_1, ..., e_n}, builds
# the linear map e_1 -> 0, e_i -> e_{i-1}
Shift := function (A)
	local B, l;
	B := CanonicalBasis (A);
	l := [Zero (A)];
	Append (l, List ([2 .. Dimension(A)], i -> B[i-1]));
	return (LeftModuleHomomorphismByImages (A, A, B, l));
end;




# A,B - algebras; returns A\otimes B
TensorProductOfTwoAlgebras := function (A, B)
	local K, scA, scB, T, list, i, j, k, l, m, n, AotimesB;

	K := LeftActingDomain(A);
	if (K <> LeftActingDomain(B)) then
		Error ("tensor factors should be defined over the same field");
	fi;

	scA := StructureConstantsTable (CanonicalBasis (A));
	scB := StructureConstantsTable (CanonicalBasis (B));

	# structure constants of the resulting algebra
	if (
		(IsCommutative     (A) and IsCommutative     (B)) or
		(IsAnticommutative (A) and IsAnticommutative (B))
		)
	then
		T := EmptySCTable (Dimension(A) * Dimension(B), Zero(K), "symmetric");  
	elif (
		(IsCommutative     (A) and IsAnticommutative (B)) or
		(IsAnticommutative (A) and IsCommutative     (B))
		)
	then		
		T := EmptySCTable (Dimension(A) * Dimension(B), Zero(K), "antisymmetric");  
	else
		T := EmptySCTable (Dimension(A) * Dimension(B), Zero(K));  
	fi;

	# assuming (a_i) is a basis of A and (b_i) is a basis B,
	# basis of A\otimes B is a set of pairs (a_i,b_j) numbered in the reverse
	# lexicographic order (so (a_i,b_j) is at i + (j-1)*Dimension(A) place)
	# TODO: take into account (anti)commutativity of A and B
	for i in [1 .. Dimension(A)] do
		for j in [1 .. Dimension(B)] do
			for k in [1 .. Dimension(A)] do
				for l in [1 .. Dimension(B)] do
					list := [];  # auxiliary list used to form structure 
								 # constants of the resulting algebra
					for n in [1 .. Length (scA[i][k][1])] do
						for m in [1 .. Length (scB[j][l][1])] do
							Add (list, scA[i][k][2][n] * scB[j][l][2][m]);
							Add (list, scA[i][k][1][n] + (scB[j][l][1][m] - 1)*Dimension(A));
						od;
					od;
					SetEntrySCTable 
						(T, i+(j-1)*Dimension(A), k+(l-1)*Dimension(A), list);
				od;
			od;
		od;
	od;
	
	AotimesB := AlgebraByStructureConstants (K, T);
	if (IsAssociative (A) and IsAssociative (B)) then
		SetIsAssociative (AotimesB, true);
	fi;

	# current Lie algebra
	if (
		(IsAssociative (A) and IsCommutative (A) and IsLieAlgebra (B)) or
		(IsAssociative (B) and IsCommutative (B) and IsLieAlgebra (A))
		) 
	then
		SetIsLieAlgebra (AotimesB, true);
	fi;

	return (AotimesB);
end;




# tensor product of several algebras as recursive tensor product of two algebras
TensorProductOfAlgebras := function (arg)
	local T, i;
	if (Length(arg) = 0) then
		# or should we return the ground field as an 1-dimensional algebra?
		Error ("no arguments supplied to TensorProductOfAlgebras()");
	fi;
	T := arg[1];
	for i in [2 .. Length(arg)] do
		T := TensorProductOfTwoAlgebras (T, arg[i]);
	od;
	return (T);
end;




# for the tensor product T = A \otimes B of two algebras A and B,
# and a linear map D on A, construct a linear map D \otimes 1 on T
MapOtimesOne := function (T, D)
	local A, basisA, basisT, dimB, images, i, j;

	basisT := CanonicalBasis (T);
	A := Source (D);
	basisA := CanonicalBasis (A);
	dimB := Dimension(T)/Dimension(A);
	if (not IsInt (dimB)) then
		Error ("dimensions do not match");
	fi;

	images := [];
	for j in [1 .. dimB] do
		for i in [1 .. Dimension(A)] do
			Add (images, LinearCombination 
				(List ([1 .. Dimension(A)], k -> basisT[k + (j-1)*Dimension(A)]),
				Coefficients (basisA, Image (D, basisA[i]))));
		od;
	od;

	return (LeftModuleHomomorphismByImages (T, T, basisT, images));
end;




# for the tensor product T = A \otimes B of two algebras A and B,
# and a linear map D on B, construct a linear map 1 \otimes D on T
OneOtimesMap := function (T, D)
	local B, basisB, basisT, dimA, images, i, j, c;

	basisT := CanonicalBasis (T);
	B := Source (D);
	basisB := CanonicalBasis (B);
	dimA := Dimension(T)/Dimension(B);
	if (not IsInt (dimA)) then
		Error ("dimensions do not match");
	fi;

	images := [];
	for j in [1 .. Dimension(B)] do
		c := Coefficients (basisB, Image (D, basisB[j]));
		for i in [1 .. dimA] do
			Add (images, LinearCombination 
				(List ([1 .. Dimension(B)], k -> basisT[i + (k-1)*dimA]), c));
				
		od;
	od;

	return (LeftModuleHomomorphismByImages (T, T, basisT, images));
end;



# end of poisson.gap