Quelle bch.g
Sprache: unbekannt
|
|
#############################################################################
#
# 7.10.2005
# Code for the computation of coefficients of the BCH-series
# in k variables. That is, we compute an expression for z in
#
# exp(y_1)* exp(y_2) *...* exp(y_k) = exp( z )
#
# First we compute z as a polynomial in the x_i. For this purpose,
# we use the paper of Reinsch "A simple expression for the terms in the
# Baker-Campbell-Hausdorff series"
#
# These terms have now to be converted in expression containing only
# Lie brackets. It is not clear to me, how you can do this for k variables,
# but this is also not really needed.
# What is is more important for us is the case
#
# Long Commutator in exp( y_1 ),exp(y_2) = exp( z )
#
# For this we replace the y_i by y_1,-y_1,y_2,-y_2.
#
# The resulting polynomial in y_1,y_2 can be transformed into an
# expression using Lie brackets in y_1, y_2 with the trick of Vaughan-Lee
#
# The aim of all this, is the computation of structure constants and
# an efficient log and exp. We use
#
# [logx,logy]= log( [x,y] ) + higher commutators in x,y (*)
#
# Maybe it is possible to compute a series for this, i.e. we would
# get [logx,logy] = Q-linear combination of logs of group commutators.
# Approach: Compute this expression for a group of fixed nilpotency
# class. Then we can probably show that a group of nilpotency class
# c+1 contains the old series plus some new terms.
# We can compute this by replacing in (*) the higher commutators in x,y
# chi(x,y)_L by chi(x,y)_G + higher commutators in x,y, which we
# compute with the extended BCH-series.
#
# Maybe we can also find a formula for log applied to
# g_1^e_1 * ... * g_k^e_k where k varies between 1 and 10.
# What we need for this is exactly the Lie bracket expression for
# z in exp(z) = exp(y_1) * ... exp( y_n )
#
# TODO:
# - BCH series for any commutator in x,y
# (also possible with Reinsch by using F,F^-1,G,G^-1)
# - BCH in k variable can be used for efficient log.
# Can we deal with elements of the factors of the upper central
# series simultaneously ?
# If you have a group of class 7 with 60 generators, it would
# be an overkill to use the BCH formula in 60 Variables!
# - Should we use basic commutators for all series ?
#############################################################################
#
# IN: v ..... upper or lower trinagular matrix with 0s on the diagonal
#
# OUT: exp(x)
#
BCH_Exponential := function( A, v )
local n,exp, vPow, fac, i;
n := Length( v[1] );
exp := One( A )* IdentityMat(n) + v;
vPow := v;
fac := 1;
for i in [2..n-1] do
vPow := vPow*v;
fac := fac*i;
exp := exp + (fac^-1)*vPow;
od;
return exp;
end;
#############################################################################
#
# IN: x ..... upper or lower trinagular matrix with ones on the diagonal
#
# OUT: log(x)
#
BCH_Logarithm := function( x )
local n,mat,matPow,log,i;
n := Length( x[1] );
mat := x - x^0;
matPow := mat;
log := mat;
for i in [2..n-1] do
matPow := matPow*mat;
log := log + (-1)^(i-1)*(i)^-1*matPow;
od;
return log;
end;
#############################################################################
#
# matrix F as in paper of Reinsch
#
BCH_Compute_F := function( n )
local F,i,j;
F := IdentityMat( n+1 );
for i in [1..n+1] do
for j in [i+1..n+1] do
F[i][j] := 1/Factorial( j-i );
od;
od;
return F;
end;
# compute F^pow where pow = +- 1
BCH_Compute_F_Power := function( n, pow )
local mat,i,F_pow;
mat := 0*IdentityMat( n+1 );
for i in [1..n] do
mat[i][i+1] := 1;
od;
F_pow := BCH_Exponential( Rationals, pow*mat );
return F_pow;
end;
BCH_Compute_G := function( n, PolRing )
local G,i,j,prod,I;
I := IndeterminatesOfPolynomialRing( PolRing );
G := IdentityMat( n+1 );
for i in [1..n+1] do
for j in [i+1..n+1] do
prod := Product( I{[i..j-1]} );
G[i][j] := 1/Factorial( j-i )* prod;
od;
od;
return G;
end;
# compute G^pow where pow = +- 1
BCH_Compute_G_Power := function( n, PolRing, pow )
local mat,i,G_pow,I;
mat := 0*IdentityMat( n+1 );
I := IndeterminatesOfPolynomialRing( PolRing );
for i in [1..n] do
mat[i][i+1] := I[i];
od;
G_pow := BCH_Exponential( Rationals, pow*mat );
return G_pow;
end;
BCH_Compute_H := function( n, PolRing )
local H,i,j,prod,I;
I := IndeterminatesOfPolynomialRing( PolRing );
H := IdentityMat( n+1 );
for i in [1..n+1] do
for j in [i+1..n+1] do
prod := Product( I{[(i+n)..(j+n-1)]} );
H[i][j] := 1/Factorial( j-i )* prod;
od;
od;
return H;
end;
#############################################################################
#
# IN: mat ............ unipotent matrix
# pow ............ rational number
#
# OUT: mat^pow
BCH_RationalPowerOfUnipotentMat := function( mat, pow )
local log;
log := BCH_Logarithm( mat );
return BCH_Exponential( Rationals, pow*log );
end;
#############################################################################
#
# Extension of the code of Reinsch. We do computations with k variables
#
# IN:
# n ........... n-th term of the extended BCH series should be computed
# k ........... k variables, i.e. exp(y_1)*...*exp(y_k)= exp(z)
# A ........... Free associative algebra with one containing the y_i
#
# For simplification of the programming we replace the matrix F, which
# originally does not contain varialbles with one similar to G,
# i.e. we get a list of k matrices called H, whose i-th matrix
# contains only y_i
#
#A := FreeAssociativeAlgebraWithOne( Rationals, k, "x" );
#
BCH_Compute_HMatrices := function( n, k, A )
local gens,i,j,mat,H;
gens := GeneratorsOfAlgebraWithOne( A );
H := [];
for i in [1..k] do
mat := Zero( A ) * IdentityMat( n+1 );
for j in [1..n] do
mat[j][j+1] := gens[i];
od;
H[i] := BCH_Exponential( A, mat );
od;
return H;
end;
# typical input string = "1/180*x_1*x_2*x_5"
BCH_SigmaExpression2XY := function( n, string, A )
local ll,k,gens,monomial,start,ll2,i,rat;
ll := SplitString( string, "*" );
k := Length( ll );
#check if start with rational number
rat := 1;
start := 1;
if ll[1][1] <> 'x' then;
rat := Rat( ll[1] );
start := 2;
fi;
ll2 := List( [start..k], x-> Int( ll[x]{[3]}) );
# compute monomial
gens := GeneratorsOfAlgebra( A );
monomial := gens[1];
for i in [1..n] do
if i in ll2 then
monomial := monomial*gens[3];
else
monomial := monomial*gens[2];
fi;
od;
return rat*monomial;
end;
#############################################################################
#
# Operator T corresponding to Reinsch
# IN
# string .......... polynomial in sigmaexpressions
# typical input string = "1/180*x_1*x_2*x_5"
# OUT
# result of operator T.
#
# Example:
# gap> BCH_SigmaExpression2List( 7, "1/180*x_1*x_2*x_5" );
# [ 1/180, [ 2, 2, 1, 1, 2, 1, 1 ] ]
#
BCH_SigmaExpression2List := function( n, string )
local ll,k,start,ll2,i,rat,result;
ll := SplitString( string, "*" );
k := Length( ll );
#check if start with rational number or minus sign
rat := 1;
start := 1;
if ll[1][1] <> 'x' then
if ll[1][1] = '-' and ll[1][2] = 'x' then
rat := -1;
ll[1] := ll[1]{[2..Length(ll[1])]};
else
rat := Rat( ll[1] );
start := 2;
fi;
fi;
# extract indices
ll2 := List( [start..k],
x-> Int( ll[x]{[3..Length(ll[x])]}) );
# compute monomial
result := [];
for i in [1..n] do
if i in ll2 then
# the 1 stands for y
Add( result, 2);
else
Add( result, 1);
fi;
od;
return [rat,result];
#return [result, rat ];
end;
#############################################################################
#
# IN:
# string .............. represents a MONOMIAL of a Free Algebra with One
# for example "(4/5)*y.2*y.1"
# OUT
# A list describing this element, in the example [ 4/5, [ 2, 1 ] ]
#
BCH_Monomial2List := function( string )
local ll,b,a,rat,inds,index,i,power,j,ss;
# get rid of new lines commands
ll := ReplacedString( string, "\n", "" );
ll := SplitString( ll, "*" );
b := Length( ll );
# get coefficient
a := Length( ll[1] );
rat := Rat( ll[1]{[2..a-1]} );
# get list which describes monomial
inds := [];
for i in [2..b] do
# deal with powers
power := 1;
ss := SplitString( ll[i], "^" );
if Length( ss ) > 1 then
power := Int( ss[2] );
fi;
a := Length( ss[1] );
index := Int( ss[1]{[3..a]} );
for j in [1..power] do
Add( inds, index );
od;
od;
return [ rat, inds ];
end;
#############################################################################
#
# A := FreeAlgebraWithOne( Rationals, 2 );
# typical input string = "-1/720*x_1*x_2+1/180*x_3*x_5";
# Conversion corresponding to Reinsch
# Output lies in A
#
BCH_SumOfSigmaExpresions2XY := function( n, string, A )
local ll1,ll2,ll3,k, gens,pol,i;
ll1 := ReplacedString( string, "-", "+-" );
ll2 := SplitString( ll1, "+" );
ll3 := Filtered( ll2, x-> x <> "" );
k := Length( ll3 );
gens := GeneratorsOfAlgebra( A );
pol := gens[1]*0;
for i in [1..k] do
pol := pol + BCH_SigmaExpression2XY( n, ll3[i], A );
od;
return pol;
end;
#############################################################################
#
# typical input string = "-1/720*x_1*x_2+1/180*x_3*x_5";
# Conversion corresponding to Reinsch
#
BCH_SumOfSigmaExpresions2List := function( n, string )
local ll1,ll2,ll3,k,i,result;
ll1 := ReplacedString( string, "-", "+-" );
ll2 := SplitString( ll1, "+" );
ll3 := Filtered( ll2, x-> x <> "" );
k := Length( ll3 );
result := [];
for i in [1..k] do
Add( result , BCH_SigmaExpression2List( n, ll3[i] ));
od;
return result;
end;
#############################################################################
#
# IN:
# string .............. represents a POLYNOMIAL of a Free Algebra with One
# for example "(1/2)*y.1*y.2+(-1/2)*y.2*y.1"
# OUT
# A list describing this element, in the example
# [ [ 1/2, [ 1, 2 ] ], [ -1/2, [ 2, 1 ] ] ]
#
BCH_Polynomial2List := function( string )
local ll,y,result,i;
ll := SplitString( string, "+" );
y := Length( ll );
result := [];
for i in [1..y] do
Add( result , BCH_Monomial2List( ll[i] ));
od;
return result;
end;
BCH_WeightInY := function( list )
local weight,i;
weight := 0;
for i in list do
if i = 2 then
weight := weight + 1;
fi;
od;
return weight;
end;
#############################################################################
#
#
# IN:
# ll ............ a list (for example [ -1/720, [ 1, 2 ] ] ) which
# specifies a monomial in x,y over Q
# Convention: 1 stands for x
# 2 stands for y
# OUT:
# a list which specifies an elment of the Lie algebra, written
# with Lie brackets, which is the image under the map of Vaughan-Lee
#
BCH_MonomialList2LieBracketList := function( ll )
local k,result,i,l,sum,weight;
k := Length( ll );
result := [];
for i in [1..k] do
# take only terms starting with yx
if ll[i][2][1]=2 and ll[i][2][2]=1 then
l := StructuralCopy( ll[i] );
# weight in y
weight := BCH_WeightInY( l[2] );
#sum := Sum( l[2] );
l[1] := l[1]/weight;
Add( result, l );
fi;
od;
return result;
end;
#############################################################################
#
#
# IN:
# ll ............ a list (for example [ -1/720, [ 1, 2 ] ] ) which
# specifies a monomial in y_1,...,y_k over Q
# Convention:
# 1 stands for y_1
# 2 stands for y_2
# and so on
# OUT:
# a list which specifies an elment of the Lie algebra, written
# with Lie brackets, which is the image under the map of Jacobson:
# d(m) = 1/i [m]
# where "i" is the weigth of the monomial.
#
# TODO Is this function really correct for more than 2 variables ?
#
BCH_MonimalInKVariablesList2LieBracketList := function( ll )
local k,result,i,l,weight;
k := Length( ll );
result := [];
for i in [1..k] do
l := StructuralCopy( ll[i] );
weight := Length( l[2] );
l[1] := l[1]/weight;
Add( result, l );
od;
return result;
end;
#############################################################################
#
# Calculation of the coefficients of the BCH series of weigt n.
#
# Eample:
#
# gap> BCH_ComputeCoefficientsAndBrackets( 2 );
# [ [ -1/2, [ 1, 0 ] ] ]
# gap> BCH_ComputeCoefficientsAndBrackets( 3 );
# [ [ -1/12, [ 1, 0, 1 ] ], [ 1/12, [ 1, 0, 0 ] ] ]
#
BCH_ComputeCoefficientsAndBrackets := function( n )
local P,F,G,elm,string,ll,lie ;
P := PolynomialRing( Rationals, n );
F := BCH_Compute_F( n );;
G := BCH_Compute_G( n, P );;
elm := BCH_Logarithm( F*G )[1][n+1];
string := String( elm );
ll := BCH_SumOfSigmaExpresions2List( n, string );
lie := BCH_MonomialList2LieBracketList( ll );
return lie;
end;
# compute bch series up to weight c
BCH_Series := function( c )
local sers, terms,i;
sers := [];
for i in [2..c] do
terms := BCH_ComputeCoefficientsAndBrackets( i );
sers := Concatenation( sers, terms );
od;
return sers;
end;
# compute bch series up to weight c
BCH_SeriesOrdered := function( c )
local sers, terms,i;
sers := [];
for i in [2..c] do
terms := BCH_ComputeCoefficientsAndBrackets( i );
Add( sers, terms );
od;
return sers;
end;
BCH_ComputeMonomials := function( n )
local P,F,G,elm,string,A,ll,lie ;
P := PolynomialRing( Rationals, n );
F := BCH_Compute_F( n );;
G := BCH_Compute_G( n, P );;
elm := BCH_Logarithm( F*G )[1][n+1];
string := String( elm );
#A := FreeAlgebraWithOne( Rationals, 2 );
ll := BCH_SumOfSigmaExpresions2List( n, string );
return ll;
end;
#############################################################################
#
# IN:
# n ......... weight, must be bigger than 1.
# word ...... list specifying word in e^y_1,e^y_2
# For example [-1,-2,1,2] means y_1^-1 * y_2^-1* y_1 * y_2
#
# Compute z where e^z = word( e^x, e^y ).
# Note that word can be any product of e^x,e^y,e^-x,e^-y
#
BCH_ComputeCoefficientsAndBracketWord := function( n, word )
local P,F,G,elm,string,ll,lie,mats,i,prod,l ;
P := PolynomialRing( Rationals, n );
mats := [];
# the following could be done better for longer words
# by saving the F,F^-1,G,G^-1
l := Length( word );
for i in [1..l] do
if AbsoluteValue( word[i] ) = 1 then
Add( mats, BCH_Compute_F_Power( n, SignInt( word[i] ) ));
elif AbsoluteValue( word[i] ) = 2 then
Add( mats, BCH_Compute_G_Power( n, P, SignInt( word[i] ) ));
else
Error( "Wrong Input\n" );
fi;
od;
prod := Product( mats );
elm := BCH_Logarithm( prod )[1][n+1];
string := String( elm );
ll := BCH_SumOfSigmaExpresions2List( n, string );
lie := BCH_MonomialList2LieBracketList( ll );
return lie;
end;
BCH_InverseWord := function( word )
local res;
res := Reversed( word );
res := -1*res;
return res;
end;
BCH_Comm2Word := function( comm )
local l,res,word,wordInv;
l := Length( comm );
if l = 2 then
res := [ -comm[1],-comm[2],comm[1],comm[2] ];
else
word := BCH_Comm2Word( comm{[1..l-1]} );
wordInv := BCH_InverseWord( word );
res := Concatenation( wordInv, [-comm[l]], word, [comm[l]] );
fi;
return res;
end;
#############################################################################
#
# IN:
# n ......... weight, must be bigger than 1.
# kappa ..... list specifyinga commutator in e^y_1,e^y_2
# For example [1,2] means [y_1,y_2]
#
# Compute z where e^z = kappa( e^x, e^y ).
#
BCH_ComputeCoefficientsAndBracketsByCommutator := function( n, kappa )
local word;
# check if
if Length( kappa ) = n then
return [ [1,kappa] ];
else
word := BCH_Comm2Word( kappa );
return BCH_ComputeCoefficientsAndBracketWord( n, word );
fi;
end;
# BCH in k variables
BCH_ComputeMonomialsExtendedBCHSeries := function( n, k )
local A,H,mat,elm,string,ll;
A := FreeAssociativeAlgebraWithOne( Rationals, k, "y" );
H := BCH_Compute_HMatrices( n, k, A );
mat := Product( H );
elm := BCH_Logarithm( mat )[1][n+1];
string := StringPrint( elm );
ll := BCH_Polynomial2List( string );
return ll;
end;
# BCH in k variables
# e^z = e^y_1* ... * e^y_k
#
BCH_ComputeCoefficientsAndBracketsExtendedBCH := function( n, k )
local A,H,mat,elm,string,ll,lie;
A := FreeAssociativeAlgebraWithOne( Rationals, k, "y" );
H := BCH_Compute_HMatrices( n, k, A );
mat := Product( H );
elm := BCH_Logarithm( mat )[1][n+1];
string := StringPrint( elm );
ll := BCH_Polynomial2List( string );
lie := BCH_MonimalInKVariablesList2LieBracketList( ll );
return lie;
end;
#############################################################################
#
# Compute all terms of the BCH series up to a given weight wCom
# of all commutators up to weight wSers.
# That is, we compute for all group commutators
# kappa of weight up to w2, all terms of weight up to w1 in the expression
# z where
# kappa( e^x, e^y ) = e^z
#
#
BCH_ComputeCommutatorSeries := function( wSers, wCom )
local coms,i,a,lie,j,n,com,series,terms;
if wSers < wCom then
Error( "The BCH series of a commutator of Length wCom starts with\n",
"Lie brackets of Length wCom.\n",
"Therefore wSers should be at least as big as wCom." );
fi;
# produce a list coms, containing all possible commutators up
# to weight wCom
coms := [];
coms[1] := -1;
coms[2] := [ [1,2], [2,1] ];
for i in [3..wCom] do
coms[i] := [];
for a in coms[i-1] do
Add( coms[i], Concatenation( a, [1] ) );
Add( coms[i], Concatenation( a, [2] ) );
od;
od;
# for each commutator compute its BCH series up to weight wSers
lie := [];
lie[1] := -1;
for i in [2..wCom] do
lie[i] := [];
for j in [1..Length(coms[i])] do
com := coms[i][j];
series := [];
# the BCH of com starts with weight length of com
for n in [Length(com)..wSers] do
terms := BCH_ComputeCoefficientsAndBracketsByCommutator(
n, com );
series := Concatenation( series, terms );
od;
Add( lie[i], series );
od;
od;
return rec( coms := coms, lie := lie );
end;
#############################################################################
#
#
# kappa( logx, logy ) = log( kappa(x,y) ) + tail
# factor*kappa( logx, logy ) = factor*log( kappa(x,y) ) + factor* tail
#
BCH_GetTail := function( kappa, sersOfComs, factor )
local w, pos,sers,i,l;
w := Length( kappa );
pos := Position( sersOfComs.coms[w], kappa );
sers := StructuralCopy( sersOfComs.lie[w][pos] );
l := Length( sers );
# get rid of starting term kappa( logx, logy )
sers := sers{[2..l]};
# change sign because we have to bring to the other side
# and multiply it with factor
for i in [1..l-1] do
sers[i][1] := -factor*sers[i][1];
od;
return sers;
end;
# add lie brackets that are similar.
# for example [ [1/4, [ 2, 1, 2, 1 ] ], [1/4, [ 2, 1, 2, 1 ] ]]
# becomes [ [1/2, [ 2, 1, 2, 1 ] ] ]
BCH_SimplifySers := function( sers )
local res,l,term,contained,term2,pos,i,res_without_zeros;
res := [];
l := Length( sers );
if l = 0 then
return [];
fi;
Add( res, sers[1] );
for i in [2..l] do
term := sers[i];
# test if term[2] is already contained in res
contained := false;
for pos in [1..Length(res)] do
term2 := res[pos];
if term2[2] = term[2] then
contained := true;
break;
fi;
od;
if contained then
# summarize
res[pos][1] := res[pos][1] + term[1];
else
# else add term
Add( res, term );
fi;
od;
# delete entries with zero coeff
res_without_zeros := [];
for pos in [1..Length( res )] do
# if zero coeff then add
if res[pos][1] <> 0 then
Add( res_without_zeros, res[pos] );
fi;
od;
return res_without_zeros;
end;
#############################################################################
#
# Aim:
#[logx,logy] = Q-linear combination of logs of group commutators.
#
#
# c .............. class of nilpotent groups for which the result
# of the computation should be used.
# That is, c is the maximal length of a group commutator
# which is non-trivial.
#
# Approach: Compute this expression for a group of fixed nilpotency
# class. Then we can probably show that a group of nilpotency class
# c+1 contains the old series plus some new terms.
# We can compute this by replacing in (*) the higher commutators in x,y
# chi(x,y)_L by chi(x,y)_G + higher commutators in x,y, which we
# compute with the extended BCH-series.
#
#
BCH_LieBracketInTermsOfLogs := function( c )
local sersOfComs,logs,term,lieTerms,tail,i,logs_ordered,le;
# get bch series of all commutators up to weight c
sersOfComs := BCH_ComputeCommutatorSeries( c, c );
# we know the first term of the logs expressions
logs := [ [1,[1,2]] ];
# get remaining terms in logx, logy
lieTerms := BCH_GetTail( [1,2], sersOfComs, 1 );
# replace iteratively all remaining Lie bracktes
i := 1;
while i <= Length( lieTerms ) do
term := lieTerms[i];
# replace it by factor*log( kappa(x,y) ) + factor*tail
Add( logs, term );
tail := BCH_GetTail( term[2], sersOfComs, term[1] );
Append( lieTerms, tail );
i := i+1;
od;
logs := BCH_SimplifySers( logs );
# order result corresponding to weight
logs_ordered := [];
for i in [2..c] do
logs_ordered[i-1] := [];
od;
for term in logs do
le := Length( term[2] );
Add( logs_ordered[le-1], term );
od;
return logs_ordered;
end;
BCH_TGroupRec := function( N )
local sers,i, NN,l,leFactors,s,indices,class,
wei, weights,a, largestAbelian;
sers := UpperCentralSeries( N );
NN := PcpGroupBySeries( sers );
#basisNN := GeneratorsOfGroup( NN );
#iso := NN!.bijection;
#basisN := GeneratorsOfGroup( Image( iso ) );
# get indices of generators of the factors of the upper central series
l := Length( sers );
leFactors := [];
for i in [1..l-1] do
Add( leFactors, HirschLength( sers[i]) -
HirschLength( sers[i+1]));
od;
s := 1;
indices := [];
for i in [1..l-1] do
Add( indices, [s..s+leFactors[i]-1] );
s := s + leFactors[i];
od;
# nilpotency class of group
class := Length( leFactors );
# weights of generators
wei := 1;
weights := [];
for a in leFactors do
for i in [1..a] do
Add( weights, wei );
od;
wei := wei + 1;
od;
# get larges abelian group in the upper central series
for i in [1..Length(sers)] do
if IsAbelian( sers[i] ) then
largestAbelian := i;
break;
fi;
od;
return rec( N := N, #basisN := basisN,
NN := NN,
sers := sers, indices := indices,
class := class, weights := weights,
largestAbelian := largestAbelian );
end;
BCH_SetUpLieAlgebraRecordByMalcevbasis := function( recTGroup )
local NN,hl,T,L;
# get dimension of algebra
NN := recTGroup.NN;
hl := HirschLength( NN );
# get prototype for structure constant table and Lie algebra
T:= EmptySCTable( hl, 0, "antisymmetric" );
L:= LieAlgebraByStructureConstants( Rationals, T );
# set the structure constants we already know,
# i.e. lie brackets involving log g_i where g_i in Z(NN)
# not necessary, because empty entries are interpreted
# as commuting elements
return rec( L := L, recTGroup := recTGroup, scTable := T );
end;
BCH_FullBCHInformation := function( w )
local bchSers, bchLBITOL;
bchSers := BCH_SeriesOrdered( w );
bchLBITOL := BCH_LieBracketInTermsOfLogs( w );
return rec( bchSers := bchSers,
bchLBITOL := bchLBITOL );
end;
# x corresponds to 1
# y corresponds to 2
BCH_EvaluateLieBracket := function( x, y, com, info )
local r,l,tmp,i;
tmp := [x,y];
r := tmp[com[1]];
l := Length( com );
if info = "lieAlgebra" then
for i in [2..l] do
r := r * tmp[com[i]];
od;
elif info = "matrix" then
for i in [2..l] do
r := LieBracket( r, tmp[com[i]] );
od;
else
Error( "wront input \n " );
fi;
return r;
end;
# g corresponds to 1
# h corresponds to 2
BCH_EvaluateGroupCommutator := function( g, h, com )
local tmp,r,l,i;
tmp := [g,h];
r := tmp[com[1]];
l := Length( com );
for i in [2..l] do
r := Comm( r, tmp[com[i]] );
od;
return r;
end;
BCH_WeightOfCommutator := function( com, wx, wy )
local weight,a;
weight := 0;
for a in com do
if a = 1 then
weight := weight + wx;
elif a = 2 then
weight := weight + wy;
else
Error( "Wrong input\n" );
fi;
od;
return weight;
end;
BCH_CheckWeightOfCommutator := function( com, wx, wy, class )
local w;
w := BCH_WeightOfCommutator( com, wx, wy );
if w > class then
return false;
else
return true;
fi;
end;
# Very SIMPLE Implementation !!
#
# compute x*y, where "*" is the BCH operation
# wx ..... weight of x
# wy ..... weight of y
# class ...... nilpotency class of the lie algebra,
# so brackets of Length class + 1 are always 0.
#
# info ..... strings which determines how lie bracket is evaluated
#
BCH_Star_Simple := function( recBCH, x, y, wx, wy, class, info )
local i,r,bchSers,com,a,term,max,min,bound;
bchSers := recBCH.bchSers;
# start with terms which are not given by Lie brackets
r := x + y;
# trivial check
if x = 0*x or y = 0*y then
return r;
fi;
# compute upper bound for the Length of commutators, which
# can be involved
max := Maximum( wx,wy );
min := Minimum( wx,wy );
# max + min* (bound-1 ) <= class
bound := Int( (class-max)/min + 1 );
# up to bound compute the commutators and add them.
# Note that the list contains commutators of length i at position i-1.
for i in [1..bound-1] do
for term in bchSers[i] do
com := term[2];
# check if weight of commutator is not to big
if BCH_CheckWeightOfCommutator( com, wx, wy, class ) then
a := BCH_EvaluateLieBracket( x, y, com, info );
r := r + term[1]*a;
fi;
od;
od;
return r;
end;
BCH_WeightOfLieAlgElm := function( recLieAlg, elm )
local basisL,coeff,i,e;
basisL := Basis( recLieAlg.L );
coeff := Coefficients( basisL, elm );
for i in [1..Length(coeff)] do
e := coeff[i];
if e <> 0 then
return recLieAlg.recTGroup.weights[i];
fi;
od;
return recLieAlg.recTGroup.class;
end;
BCH_AbstractLog_Simple_ByExponent := function( recLieAlg, recBCH, exp )
local gensL,r,i,e,x,y,wx,wy,class;
# setup
gensL := GeneratorsOfAlgebra( recLieAlg.L );
class := recLieAlg.recTGroup.class;
r := Zero( recLieAlg.L );
# go backwards through exponents.
# Which direction is more efficient ?
for i in Reversed([1..Length( exp )] ) do
e := exp[i];
if e <> 0 then
x := e*gensL[i];
y := r;
# if y is the identity then x*y = x
if y = Zero( recLieAlg.L ) then
r := x;
else
# compute weights of x,y
wx := recLieAlg.recTGroup.weights[i];
wy := BCH_WeightOfLieAlgElm( recLieAlg, y );
r := BCH_Star_Simple( recBCH,x,y,wx,wy,class, "lieAlgebra" );
fi;
fi;
od;
return r;
end;
BCH_AbstractLogCoeff_Simple_ByExponent := function( recLieAlg, recBCH, exp )
local r;
r := BCH_AbstractLog_Simple_ByExponent( recLieAlg, recBCH, exp );
return Coefficients( Basis( recLieAlg.L ), r);
end;
BCH_AbstractLog_Simple_ByElm := function( recLieAlg, recBCH, g )
local exp,hl,l,expNN;
# get exponets with respect to gens of T group NN
# (g might an elment of PcpGroup which contains NN as a normal sugroup,
# then Exponents(g) would be longer hl(NN))
exp := Exponents( g );
hl := HirschLength( recLieAlg.recTGroup.NN );
l := Length( exp );
expNN := exp{[l-hl+1..l]};
return BCH_AbstractLog_Simple_ByExponent( recLieAlg, recBCH, expNN );
end;
BCH_AbstractLogCoeff_Simple_ByElm := function( recLieAlg, recBCH, g )
local r;
r := BCH_AbstractLog_Simple_ByElm( recLieAlg, recBCH, g );
return Coefficients( Basis( recLieAlg.L ), r);
end;
BCH_EvaluateLieBracketsInTermsOfLogarithmsSers
:= function( recBCH, recLieAlg, g,h,wg,wh )
local bchLBITOL,r,max,min,bound,class,i,term,com,a,log_a;
bchLBITOL := recBCH.bchLBITOL;
r := Zero( recLieAlg.L );
# compute upper bound for the Length of commutators, which
# can be involved
class := recLieAlg.recTGroup.class;
max := Maximum( wg,wh );
min := Minimum( wg,wh );
# max + min* (bound-1 ) <= class
bound := Int( (class-max)/min + 1 );
# up to bound compute the commutators and add them.
# Note that the list contains commutators of length i at position i-1.
for i in [1..bound-1] do
for term in bchLBITOL[i] do
com := term[2];
# check if weight of commutator is not to big
if BCH_CheckWeightOfCommutator( com, wg, wh, class ) then
# evaluate commutator in the group
a := BCH_EvaluateGroupCommutator( g, h, com );
# map to the Lie algebra
log_a := BCH_AbstractLog_Simple_ByElm( recLieAlg, recBCH, a );
r := r + term[1]*log_a;
fi;
od;
od;
return r;
end;
# output is a list, as required for SetEntrySCTable
BCH_LieAlgElm2CoeffGenList := function( L, x )
local basis, coeff,i,ll;
basis := Basis( L );
coeff := Coefficients( basis, x );
ll := [];
for i in [1..Length(coeff)] do
if coeff[i] <> 0 then
Append( ll, [coeff[i],i] );
fi;
od;
return ll;
end;
#recBCH := BCH_FullBCHInformation( 5 );
#N := PcpGroupByMatGroup( PolExamples(4 ) );
#recTGroup := BCH_TGroupRec( N );
#recLieAlg := BCH_SetUpLieAlgebraRecordByMalcevbasis( recTGroup );
# F := FreeGroup( 3 );
# N := NilpotentQuotient( F, 5 );
# recTGroup := BCH_TGroupRec( N );
# recLieAlg := BCH_SetUpLieAlgebraRecordByMalcevbasis( recTGroup );
BCH_ComputeStructureConstants := function( args )
local factors, index_x, index_y, indices,indicesNoCenter,l,g,h,wg,wh,
lie_elm,ll,gens,T,out,recLieAlg,recBCH;
# setup
recLieAlg := args[2];
recBCH := args[1];
indices := recLieAlg.recTGroup.indices;
l := Length( indices );
indicesNoCenter := indices{[1..l-1]};
gens := GeneratorsOfGroup( recLieAlg.recTGroup.NN );
T := recLieAlg.scTable;
# go through blocks of generators backwards,
# where a block corresponds to the generators of the factors
# of the upper central series
for factors in Reversed( indicesNoCenter ) do
for index_y in Reversed( factors ) do
for index_x in [1..index_y-1] do
g := gens[index_x];
h := gens[index_y];
wg := recLieAlg.recTGroup.weights[index_x];
wh := recLieAlg.recTGroup.weights[index_y];
# compute [log(g),log(h)]
lie_elm := BCH_EvaluateLieBracketsInTermsOfLogarithmsSers(
recBCH, recLieAlg, g,h, wg, wh );
# set entry in structure constant table
ll := BCH_LieAlgElm2CoeffGenList( recLieAlg.L, lie_elm );
if Length( ll ) > 0 then
SetEntrySCTable( T, index_x, index_y, ll );
fi;
od;
od;
# update Lie algebra with new structure constant table
recLieAlg.L:= LieAlgebraByStructureConstants( Rationals, T );
od;
return 0;
end;
BCH_Abstract_Exponential_ByElm := function( recBCH, recLieAlg, x )
local indices,basis,class,tail,coeffs,largestAbelian,exp_x,i,factor,
divider,w_divider,w_tail,l,exp_x_2ndPart,j;
indices := recLieAlg.recTGroup.indices;
basis := Basis( recLieAlg.L );
class := recLieAlg.recTGroup.class;
# divide some exponents untill the remaining element lies in an
# abelian group.
tail := x;
coeffs := Coefficients( basis, tail );
largestAbelian := recLieAlg.recTGroup.largestAbelian;
exp_x := [];
for i in [1..largestAbelian-1] do
factor := indices[i];
for j in factor do
# get element to divide of
divider := -coeffs[j]*basis[j];
# save exponents of divider
Add( exp_x, coeffs[j] );
# divide off
w_divider := BCH_WeightOfLieAlgElm ( recLieAlg, divider );
w_tail := BCH_WeightOfLieAlgElm ( recLieAlg, tail );
tail := BCH_Star_Simple( recBCH, divider, tail, w_divider,
w_tail, class, "lieAlgebra" );
# set up coefficient vector
coeffs := Coefficients( basis, tail );
od;
od;
# test intermediate result
l := Length( exp_x );
if not coeffs{[1..l]} = 0 * coeffs{[1..l]} then
Error( "Failure in Abstract_Exponential \n" );
fi;
# get the remaining coefficients
exp_x_2ndPart := coeffs{[l+1..Length(coeffs)]};
return Concatenation( exp_x, exp_x_2ndPart );
end;
BCH_Abstract_Exponential_ByVector := function( recBCH, recLieAlg, vec )
local basis,x;
basis := Basis( recLieAlg.L );
x := LinearCombination( basis, vec );
return BCH_Abstract_Exponential_ByElm( recBCH, recLieAlg, x );
end;
BCH_LieAlgebraByTGroupRec := function( recBCH, recTGroup )
local recLieAlg;
recLieAlg := BCH_SetUpLieAlgebraRecordByMalcevbasis( recTGroup );
BCH_ComputeStructureConstants( [recBCH, recLieAlg] );
return recLieAlg;
end;
#############################################################################
#
# TEST Functions
#
#############################################################################
BCH_Test_LogOfExp := function( recBCH, recLieAlg, noTests )
local x,exp,x2,i;
for i in [1..noTests] do
x := Random( recLieAlg.L );
exp := BCH_Abstract_Exponential_ByElm( recBCH, recLieAlg, x );
x2 := BCH_AbstractLog_Simple_ByExponent( recLieAlg, recBCH, exp );
if not x = x2 then
Error( "Mist\n" );
fi;
od;
return 0;
end;
BCH_Test_ExpOfLog := function( recBCH, recLieAlg, noTests,range )
local i,hl,domain,exp,x,exp2;
hl := HirschLength( recLieAlg.recTGroup.NN );
domain := [-range..range];
for i in [1..noTests] do
exp := List( [1..hl], x -> Random( domain ) );
x := BCH_AbstractLog_Simple_ByExponent( recLieAlg, recBCH, exp );
exp2 := BCH_Abstract_Exponential_ByElm( recBCH, recLieAlg, x );
if not exp = exp2 then
Error( "Mist\n" );
fi;
od;
return 0;
end;
BCH_Test := function( n )
local P,F,G,elm,string,A,ll,lie,kk,kk2;
P := PolynomialRing( Rationals, n );
F := BCH_Compute_F( n );;
G := BCH_Compute_G( n, P );;
elm := BCH_Logarithm( F*G )[1][n+1];
string := String( elm );
A := FreeAlgebraWithOne( Rationals, 2 );
ll := BCH_SumOfSigmaExpresions2List( n, string );
lie := BCH_MonomialList2LieBracketList( ll );
Length( ll );
Length( lie );
kk := List( lie, x->x[2] );
kk2 := AsSet( kk );
return [ Length( lie ), Length( AsSet( kk ) )];
end;
BCH_List2Number := function( base, ll )
local n,l,numb,i;
n := Length( ll );
#throw away the first two entries 1 and 0
l := ll{[3..n]};
numb := 0;
for i in Reversed( [3..n]) do
if ll[i]<>0 then
numb := numb + base^(n-i);
fi;
od;
return numb;
end;
# produce random element of Tr_0(dim,Q)
BCH_RandomNilpotentMat := function( dim )
local range,ll,kk,g,j,k;
range := 7;
ll := [ - range .. range ];
kk := [ - range .. range ];
ll := Filtered( ll, function ( x )
return x <> 0;
end );
kk := Filtered( kk, function ( x )
return x <> 0;
end );
g := NullMat( dim, dim, Rationals );
for j in [ 1 .. dim ] do
for k in [ j .. dim ] do
if j < k then
g[j][k] := RandomList( ll ) / RandomList( kk );
else
g[j][k] := 0;
fi;
od;
od;
return g;
end;
BCH_TestBchSeriesOrdered := function( recBCH, n )
local no_tests,i,x,y,wx,wy,x_star_y,exp_x_star_y,exp_x,exp_y;
no_tests := 100;
for i in [1..no_tests] do
# produce two random matrices x,y in Tr_0(n,Q)
x := BCH_RandomNilpotentMat( n );
y := BCH_RandomNilpotentMat( n );
wx := 1;
wy := 1;
# compute exp(x*y) with BCH,
# note that we need terms of length at most n-1
x_star_y := BCH_Star_Simple ( recBCH, x, y, wx, wy, n-1, "matrix" );
exp_x_star_y := BCH_Exponential( Rationals, x_star_y );
# compute exp(x)exp(y) and compare
exp_x := BCH_Exponential( Rationals, x );
exp_y := BCH_Exponential( Rationals, y );
if not exp_x_star_y = exp_x*exp_y then
Error( "Mist \n " );
fi;
od;
return 0;
end;
# recComSers as computed by BCH_ComputeCommutatorSeries
# recComSers := BCH_ComputeCommutatorSeries( 6, 3 );
BCH_Test_ComputeCommutatorSeries := function( recComSers, n, kappa )
local no_tests,i,x,y,wx,wy,exp_x,exp_y,exp_z,r,class,
sers,com,a,term,exp_z_bch,weight, pos;
no_tests := 10;
for i in [1..no_tests] do
# produce two random matrices x,y in Tr_0(n,Q)
x := BCH_RandomNilpotentMat( n );
y := BCH_RandomNilpotentMat( n );
wx := 1;
wy := 1;
# compute kappa( exp(x),exp(y) ) normally
exp_x := BCH_Exponential( Rationals, x );
exp_y := BCH_Exponential( Rationals, y );
exp_z := BCH_EvaluateGroupCommutator( exp_x, exp_y, kappa );
# compute z where exp(z)= kappa( exp(x),exp(y) ) with extended BCH
r := NullMat( n,n );
class := n-1;
weight := Length( kappa );
pos := Position( recComSers.coms[weight], kappa );
sers := recComSers.lie[weight][pos];
for term in sers do
com := term[2];
#Print( "term ", term, "\n" );
# check if weight of commutator is not to big
if BCH_CheckWeightOfCommutator( com, wx, wy, class ) then
# evaluate commutator in the lie algebra
a := BCH_EvaluateLieBracket( x, y, com, "matrix" );
r := r + term[1]*a;
# Print( "log_a ", log_a, "\n" );
# Print( "r ", r, "\n\n" );
fi;
od;
exp_z_bch := BCH_Exponential( Rationals, r );
# compare
if not exp_z_bch = exp_z then
Error( "Mist\n" );
fi;
od;
return 0;
end;
BCH_TestSeriesLieBracketInTermsOfLogs := function( recBCH, n )
local no_tests,i,x,y,x_bracket_y,g,h,wg,wh,r,class,max,min,bound,j,
term,a,log_a,x_bracket_y_2,bchLBITOL,com;
no_tests := 10;
for j in [1..no_tests] do
# produce two random matrices x,y in Tr_0(n,Q)
x := BCH_RandomNilpotentMat( n );
y := BCH_RandomNilpotentMat( n );
# compute [x,y] in normal way
x_bracket_y := LieBracket( x, y );
# compute [x,y] with series "liebrackets in terms of logs"
g := BCH_Exponential( Rationals, x );
h := BCH_Exponential( Rationals, y );
wg := 1;
wh := 1;
bchLBITOL := recBCH.bchLBITOL;
r := NullMat( n,n );
# compute upper bound for the Length of commutators, which
# can be involved
class := n-1;
max := Maximum( wg,wh );
min := Minimum( wg,wh );
# max + min* (bound-1 ) <= class
bound := Int( (class-max)/min + 1 );
# up to bound compute the commutators and add them.
# Note that the list contains comms of length i at position i-1.
for i in [1..bound-1] do
for term in bchLBITOL[i] do
com := term[2];
#Print( "term ", term, "\n" );
# check if weight of commutator is not to big
if BCH_CheckWeightOfCommutator( com, wg, wh, class ) then
# evaluate commutator in the group
a := BCH_EvaluateGroupCommutator( g, h, com );
# map to the Lie algebra
log_a := BCH_Logarithm( a );
r := r + term[1]*log_a;
#Print( "log_a ", log_a, "\n" );
#Print( "r ", r, "\n\n" );
fi;
od;
od;
x_bracket_y_2 := r;
# compare
if not x_bracket_y = x_bracket_y_2 then
Error( "Mist\n" );
fi;
od;
return 0;
end;
[ Dauer der Verarbeitung: 0.56 Sekunden
(vorverarbeitet)
]
|
2026-04-02
|
