Quelle  dt.g   Sprache: unbekannt

#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Wolfgang Merkwitz.
##
##  Copyright of GAP belongs to its developers, whose names are too numerous
##  to list here. Please refer to the COPYRIGHT file for details.
##
##  SPDX-License-Identifier: GPL-2.0-or-later
##
##  This file initializes Deep Thought.
##
##  Deep Thought deals with trees.  A tree < tree > is a concatenation of
##  several nodes where each node is a 5-tuple of immediate integers.  If
##  < tree > is an atom it contains only one node,  thus it is itself a
##  5-tuple. If < tree > is not an atom we obtain its list representation by
##
##  < tree >  :=  topnode(<tree>) concat left(<tree>) concat right(<tree>) .
##
##  Let us denote the i-th node of <tree> by (<tree>, i)  and the tree rooted
##  at (<tree>, i) by tree(<tree>, i).  Let <a> be tree(<tree>, i)
##  The first entry of (<tree>, i) is pos(a),
##  and the second entry is num(a). The third entry of (<tree>, i) gives a
##  mark.(<tree>, i)[3] = 1  means that (<tree>, i) is marked,
##  (<tree>, i)[3] = 0 means that (<tree>, i) is not marked. The fourth entry
##  of (<tree>, i) contains the number of nodes of tree(<tree>, i).  The
##  fifth entry of (<tree>, i) finally contains a boundary for
##  pos( tree(<tree>, i) ).  (<tree>, i)[5] <= 0 means that
##  pos( tree(<tree>, i) ) is unbounded.  If tree(<tree>, i) is an atom we
##  already know that pos( tree(<tree>, i) ) is unbound.  Thus we then can
##  use the fifth component of (<tree>, i) to store the side.  In this case
##  (<tree>, i)[5] = -1 means  that tree(<tree>, i) is an atom from the
##  right hand word, and (<tree>, i)[5] = -2 means that tree(<tree>, i) is
##  an atom from the left hand word.
##
##  A second important data structure deep thought deals with is a deep
##  thought monomial. A deep thought monomial g_<tree> is a product of
##  binomial coefficients with a coefficient c. Deep thought monomials
##  are represented in this implementation by formula
##  vectors,  which are lists of integers.  The first entry of a formula
##  vector is 0,  to distinguish formula vectors from trees.  The second
##  entry is the coefficient c,  and the third and fourth entries are
##  num( left(tree) ) and num( right(tree) ).  The remaining part of the
##  formula vector is a concatenation of pairs of integers.  A pair (i, j)
##  with i > 0 represents binomial(x_i, j).  A pair (0, j) represents
##  binomial(y_gen, j) when word*gen^power is calculated.
##
##  Finally deep thought has to deal with pseudo-representatives. A
##  pseudo-representative <a> is stored in list of length 4. The first entry
##  stores left( <a> ),  the second entry contains right( <a> ),  the third
##  entry contains num( <a> ) and the last entry finally gives a boundary
##  for pos( <b> ) for all trees <b> which are represented by <a>.
##




#############################################################################
##
#F  Dt_Mkavec(<pr>)
##
##  'Dt_Mkavec' returns the avec for the pc-presentation <pr>.
##
BindGlobal( "Dt_Mkavec", function(pr)
    local  i,j,vec;

    vec := [];
    vec[Length(pr)] := 1;
    for  i in [Length(pr)-1,Length(pr)-2..1]  do
        j := Length(pr);
        while  j >= 1  do
            if  j = vec[i+1]  then
                vec[i] := j;
                j := 0;
            else
                j := j-1;
                if  j < i  and  IsBound(pr[i][j])  then
                    vec[i] := j+1;
                    j := 0;
                fi;
                if  j > i  and  IsBound(pr[j][i])  then
                    vec[i] := j+1;
                    j := 0;
                fi;
            fi;
        od;
    od;
    for  i in [1..Length(pr)]  do
        if  vec[i] < i+1  then
            vec[i] := i+1;
        fi;
    od;
    return vec;
end );



#############################################################################
##
#F  Dt_IsEqualMonomial(<vec1>, <vec2>) . . . . . .  test if <vec1> and <vec2>
##                                                represent the same monomial
##
##  'Dt_IsEqualMonomial' returns "true" if <vec1> and <vec2> represent the
##  same  monomial, and "false" otherwise.
BindGlobal( "Dt_IsEqualMonomial", function(vec1,vec2)
    local  i;

    if  Length(vec1) <> Length(vec2)  then
        return false;
    fi;
    #  Since the first four entries of a formula vector doesn't contain
    #  any information about the monomial it represents,  it suffices to
    #  compare the remaining entries.
    for  i in [5..Length(vec1)]  do
        if  not vec1[i] = vec2[i]  then
            return false;
        fi;
    od;
    return true;
end );




#############################################################################
##
#F  Dt_Sort2(<vector>) . . . . . . . . . . . . . . . .  sort a formula vector
##
##  'Dt_Sort2' sorts the pairs of integers in the formula vector <vector>
##  representing the binomial coefficients such that
##  <vector>[5] < <vector>[7] < .. < vector[m-1],  where m is the length
##  of <vector>.  This is done for a easier comparison of formula vectors.
##
BindGlobal( "Dt_Sort2", function(vector)
    local  i,list1,list2;

    list1 := vector{[5,7..Length(vector)-1]};
    list2 := vector{[6,8..Length(vector)]};
    SortParallel(list1,list2);
    for  i in [1..Length(list1)]  do
        vector[ 2*i+3 ] := list1[i];
        vector[ 2*i+4 ] := list2[i];
    od;
end );


#############################################################################
##
#F  Dt_AddVecToList(<evlist>,<evlistvec>,<formvec>,<pr>)
##
##  'Dt_AddVecToList' adds the formula vector <formvec> to the list <evlist>,
##  computes the corresponding coefficient vector and adds the latter to
##  the list <evlistvec>.
##
BindGlobal( "Dt_AddVecToList", function(evlist, evlistvec, formvec, pr)
    local    i,j,k;

    Add(evlist, formvec);
    k := [];
    for  i in [1..Length(pr)]  do
        k[i] := 0;
    od;
    j := pr[ formvec[3] ][ formvec[4] ];
    #  the coefficient that the monomial represented by <formvec> has
    #  in each polynomial f_l is obtained by multiplying <formvec>[2]
    #  with the exponent which the group generator g_l has in the
    #  in the word representing the commutator of g_(formvec[3]) and
    #  g_(formvec[4]) in the presentation <pr>.
    for  i in [3,5..Length(j)-1]  do
        k[ j[i] ] := formvec[2]*j[i+1];
    od;
    Add(evlistvec, k);
end );


#############################################################################
##
#F  Dt_add( <pol>, <pols>, <pr> )
##
##  'Dt_add' adds the deep thought monomial <pol> to the list of polynomials
##  <pols>,  such that afterwards <pols> represents a simplified polynomial.
##
BindGlobal( "Dt_add", function(pol, pols, pr)
    local  j,k,rel, pos, flag;

    # first sort the deep thought monomial <pol> to compare it with the
    # monomials contained in <pols>.
    Dt_Sort2(pol);
    # then look which component of <pols> contains <pol> in case that
    # <pol> is contained in <pols>.
    pos := DT_evaluation(pol);
    if  not IsBound( pols[pos] )  then
        # create the component <pols>[<pos>] and add <pol> to it
        pols[pos] := rec( evlist := [], evlistvec := [] );
        Dt_AddVecToList( pols[pos].evlist, pols[pos].evlistvec, pol, pr );
        return;
    fi;
    flag := 0;
    for  k in [1..Length( pols[pos].evlist ) ]  do
        # look for <pol> in <pols>[<pos>] and if <pol> is contained in
        # <pols>[<pos>] then adjust the corresponding coefficient vector.
        if  Dt_IsEqualMonomial( pol, pols[pos].evlist[k] )  then
            rel := pr[ pol[3] ][ pol[4] ];
            for  j in [3,5..Length(rel)-1]  do
                pols[pos].evlistvec[k][ rel[j] ] :=
                  pols[pos].evlistvec[k][ rel[j] ] + pol[2]*rel[j+1];
            od;
            flag := 1;
            break;
        fi;
    od;
    if  flag = 0  then
        # <pol> is not contained in <pols>[<pos>] so add it to <pols>[<pos>]
        Dt_AddVecToList(pols[pos].evlist, pols[pos].evlistvec, pol, pr);
    fi;
end );


#############################################################################
##
#F  Dt_Convert(<sortedpols>)
##
##  'Dt_Convert' converts the list of formula vectors <sortedpols>. Before
##  applying <Dt_Convert> <sortedpols> is a list of records with the
##  components <evlist> and <evlistvec> where <evlist> contains deep thought
##  monomials and <evlistvec> contains the corresponding coefficient vectors.
##  <Dt_Convert> merges the <evlist>-components of the records contained
##  in <sortedpols> into one component <evlist> and the <evlistvec>-components
##  into one component <evlistvec>.
##
BindGlobal( "Dt_Convert", function(sortedpols)
    local  k,res;

    if  Length(sortedpols) = 0  then
        return  0;
    fi;
    res := rec(evlist := [],
               evlistvec :=[]);
    for  k in sortedpols  do
        Append(res.evlist, k.evlist);
        Append(res.evlistvec, k.evlistvec);
    od;
    return res;
end );


#############################################################################
##
#F  Dt_ConvertCoeffVecs(<evlistvec>)
##
##  'Dt_ConvertCoeffVecs' converts the coefficient vectors in the list
##  <evlistvec>. Before applying <Dt_ConvertCoeffVecs>, an entry
##  <evlistvec>[i][j] = k means that the deep thought monomial <evlist>[i]
##  occurs in the polynomial f_j with coefficient k. After applying
##  <Dt_ConvertCoeffVecs> a pair [j, k] occurring in <evlistvec>[i] means that
##  <evlist>[i] occurs in f_j with coefficient k.
##
BindGlobal( "Dt_ConvertCoeffVecs", function(evlistvec)
    local i,j,res;

    for  i in [1..Length(evlistvec)]  do
        res := [];
        for  j in [1..Length(evlistvec[i])]  do
            if  evlistvec[i][j] <> 0  then
                Append(res, [j, evlistvec[i][j] ]);
            fi;
        od;
        evlistvec[i] := res;
    od;
end );



#############################################################################
##
#F  CalcOrder( <word>, <dtrws> )
##
##  CalcOrder computes the order of the word <word> in the group determined
##  by the rewriting system <dtrws>
##
DeclareGlobalName( "CalcOrder" );
BindGlobal( "CalcOrder", function(word, dtrws)
    local gcd, m;

    if  Length(word) = 0  then
        return 1;
    fi;
    if  not IsBound(dtrws![PC_EXPONENTS][ word[1] ])  then
        return 0;
    fi;
    gcd := Gcd(dtrws![PC_EXPONENTS][ word[1] ], word[2]);
    m := QuoInt( dtrws![PC_EXPONENTS][ word[1] ], gcd);
    gcd := DTPower(word, m, dtrws);
    return  m*CalcOrder(gcd, dtrws);
end );


#############################################################################
##
#F  CompleteOrdersOfRws( <dtrws> )
##
##  CompleteOrdersOfRws computes the orders of the generators of the
##  deep thought rewriting system <dtrws>
##
BindGlobal( "CompleteOrdersOfRws", function(dtrws)
    local  i,j;

    dtrws![PC_ORDERS] := [];
    for  i in [dtrws![PC_NUMBER_OF_GENERATORS],dtrws![PC_NUMBER_OF_GENERATORS]-1..1]
         do
        # Print("determining order of generator ",i,"\n");
        if  not IsBound( dtrws![PC_EXPONENTS][i] )  then
            j := 0;
        elif  not IsBound( dtrws![PC_POWERS][i] )  then
            j := dtrws![PC_EXPONENTS][i];
        else
            j := dtrws![PC_EXPONENTS][i]*CalcOrder(dtrws![PC_POWERS][i], dtrws);
        fi;
        if  j <> 0  then
            dtrws![PC_ORDERS][i] := j;
        fi;
    od;
end );


#############################################################################
##
#F  Dt_RemoveHoles( <list> )
##
##  Dt_RemoveHoles removes all empty entries from <list>
##  It is similar to Compacted(), but works in-place
##
BindGlobal( "Dt_RemoveHoles", function( list )
    local  skip, i;

    skip := 0;
    i := 1;
    while  i <= Length(list)  do
        while  not IsBound(list[i])  do
            skip := skip + 1;
            i := i+1;
        od;
        list[i-skip] := list[i];
        i := i+1;
    od;
    for  i in  [Length(list)-skip+1..Length(list)]  do
        Unbind(list[i]);
    od;
end );



#############################################################################
##
#F  ReduceCoefficientsOfRws( <dtrws> )
##
##  ReduceCoefficientsOfRws reduces all coefficients of each deep thought
##  polynomial f_l modulo the order of the l-th generator.
##
BindGlobal( "ReduceCoefficientsOfRws", function(dtrws)
    local  i,j,k, pseudoreps;

    pseudoreps := dtrws![PC_DEEP_THOUGHT_POLS];
    i := 1;
    while  IsRecord(pseudoreps[i])  do
        for  j in [1..Length(pseudoreps[i].evlistvec)]  do
            for  k in [2,4..Length(pseudoreps[i].evlistvec[j])]  do
                if  IsBound( dtrws![PC_ORDERS][ pseudoreps[i].evlistvec[j][k-1] ] )
                    and  (pseudoreps[i].evlistvec[j][k] > 0  or
                          pseudoreps[i].evlistvec[j][k] <
                          -dtrws![PC_ORDERS][ pseudoreps[i].evlistvec[j][k-1] ]/2)
                    then
                    pseudoreps[i].evlistvec[j][k] :=
                      pseudoreps[i].evlistvec[j][k] mod
                      dtrws![PC_ORDERS][ pseudoreps[i].evlistvec[j][k-1] ];
                fi;
            od;
            DTCompress( pseudoreps[i].evlistvec[j] );
            if  Length( pseudoreps[i].evlistvec[j] ) = 0  then
                Unbind( pseudoreps[i].evlistvec[j] );
                Unbind( pseudoreps[i].evlist[j] );
            fi;
        od;
        Dt_RemoveHoles( pseudoreps[i].evlistvec );
        Dt_RemoveHoles( pseudoreps[i].evlist );
        i := i+1;
    od;
end );


#############################################################################
##
##  Dt_GetMax( <tree>, <number>, <pr> )
##
##  Dt_GetMax returns the maximal value for pos(tree) if num(tree) = <number>.
##
BindGlobal( "Dt_GetMax", function(tree, number, pr)
    local rel, position;

    if  Length(tree) = 5  then
        return tree[5];
    else
        if  Length(tree) = 4  then
            if  Length(tree[1]) = 4  then
                if  Length(tree[2]) = 4  then
                    rel := pr[ tree[1][3] ][ tree[2][3] ];
                else
                    rel := pr[ tree[1][3] ][ tree[2][2] ];
                fi;
            else
                if  Length(tree[2]) = 4  then
                    rel := pr[ tree[1][2] ][ tree[2][3] ];
                else
                    rel := pr[ tree[1][2] ][ tree[2][2] ];
                fi;
            fi;
        else
            rel := pr[ tree[7] ][ tree[ 5*(tree[9]+1)+2 ] ];
        fi;
        position := Position(rel, number) + 1;
        if  rel[position] < 0  or  rel[position] > 100  then
            return 0;
        else
            return rel[position];
        fi;
    fi;
end );


#############################################################################
##
#F  Dt_GetNumRight( <tree> )
##
##  Dt_GetNumRight  returns num( right( tree ) ).
##
BindGlobal( "Dt_GetNumRight", function(tree)

    if Length(tree) <> 4  then
        return  tree[ 5*(tree[9]+1)+2 ];
    fi;
    if Length(tree[2]) <> 4  then
        return  tree[2][2];
    fi;
    return  tree[2][3];
end );



###########################################################################
##
#F  Calcrepsn(<n>, <avec>, <pr>, <max>)
##
##  'Calcrepsn' returns the polynomials f_{n1},...,f_{nm} which have to be
##  evaluated when computing word*g_n^(y_n).  Here m denotes the composition
##  length of the nilpotent group G given by the presentation <pr>.  This is
##  done  by first calculating a complete system of <n>-pseudo-representatives
##  for the presentation <pr> with boundary <max>. Then this system is used
##  to get the required polynomials
##
##  If g_n is in the center of the group determined by the presentation <pr>
##  then there don't exist any representatives except for the atoms and
##  finally 0 will be returned.
##
BindGlobal( "Calcrepsn", function(n, avec, pr, max)

    local i,j,k,l,       #  loop variables
          x,y,z,a,b,c,   #  trees
          reps,          #  list of pseudo-representatives
          pols,          #  stores the deep thought polynomials
          boundary,      #  boundary for loop
          hilf,
          pos,
          start,
          max1, max2;    #  maximal values for pos(x) and pos(y)
    reps := [];
    pols := [];

    for i  in [n..Length(pr)]  do
        #  initialize reps[i] to contain representatives for the atoms
        if  i <> n  then
            reps[i] := [ [1,i,0,1,-2] ];
        else
            reps[i] := [ [1,i,0,1,-1] ];
        fi;
    od;
    #  first compute the pseudo-representatives which are also representatives
    for  i in [n..max]  do
        if  i < avec[n]  then
            boundary := i-1;
        else
            boundary := avec[n]-1;
        fi;
        #  to get the representatives of the non-atoms loop over j and k
        #  and determine the representatives for all trees <z> with
        #  num(<z>) = i,  num( left(<z>) ) = j,  num( right(<z>) ) = k
        #  Since for all 1 <= l <= m the group generated by
        #  {g_(avec[l]),..,g_m} is in the center of the group generated
        #  by {g_l,..,g_m} it suffices to loop over all
        #  j <= min(i-1, avec[n]-1). Also it is sufficient only to loop over
        #  k while avec[k] is bigger than j.
        for j in [n+1..boundary] do
            k := n;
            while  k <= j-1  and  avec[k] > j  do
                if  IsBound(pr[j][k])  and  pr[j][k][3] = i  then
                    if  k = n  then
                        start := 1;
                    else
                        start := 2;
                    fi;
                    for x in [start..Length(reps[j])] do
                        for y in reps[k] do
                            if Length(reps[j][x]) = 5
                               or  k >= reps[j][x][ 5*(reps[j][x][9]+1)+2 ]
                               then
                                max1 := Dt_GetMax(reps[j][x], j, pr);
                                max2 := Dt_GetMax(y, k, pr);
                                z := [1,i,0, reps[j][x][4]+y[4]+1, 0];
                                Append(z,reps[j][x]);
                                Append(z,y);
                                z[7] := j;
                                z[10] := max1;
                                z[ 5*(z[9]+1)+2 ] := k;
                                z[ 5*(z[9]+1)+5 ] := max2;
                                UnmarkTree(z);
                                #  now get all representatives <z'> with
                                #  left(<z'>) = left(<z>)  ( = <x> ) and
                                #  right(<z'>)=right(<z>) ( = <y> ) and
                                #  num(<z'>) = o where o is an integer
                                #  contained in pr[j][k].
                                FindNewReps(z, reps, pr, avec[n]-1);
                            fi;
                        od;
                    od;
                fi;
                k := k+1;
            od;
        od;
    od;
    #  now get the "real" pseudo-representatives
    for  i in [max+1..Length(pr)]  do
        if  i < avec[n]  then
            boundary := i-1;
        else
            boundary := avec[n]-1;
        fi;
        for j in [n+1..boundary] do
            k := n;
            while  k <= j-1  and  avec[k] > j  do
                if  IsBound(pr[j][k])  and  pr[j][k][3] = i  then
                    if  k = n  then
                        start := 1;
                    else
                        start := 2;
                    fi;
                    for x in [start..Length(reps[j])]  do
                        for y in reps[k]  do
                            if Length(reps[j][x]) = 5
                               or  k >= Dt_GetNumRight(reps[j][x])  then
                                # since reps[j] and reps[k] may contain
                                # pseudo-representatives which are trees
                                # as well as "real" pseudo-representatives
                                # it is necessary to take several cases into
                                # consideration.
                                max1 := Dt_GetMax(reps[j][x], j, pr);
                                max2 := Dt_GetMax(y, k, pr);
                                if  Length(reps[j][x]) <> 4  then
                                    if  reps[j][x][2] <> j  then
                                        # we have to ensure that
                                        # num( <reps>[j][x] ) = j when we
                                        # construct a new pseudo-representative
                                        # out of it.
                                        a := ShallowCopy(reps[j][x]);
                                        a[2] := j;
                                    else
                                        a := reps[j][x];
                                    fi;
                                    a[5] := max1;
                                else
                                    if  reps[j][x][3] <> j  then
                                        # we have to ensure that
                                        # num( <reps>[j][x] ) = j when we
                                        # construct a new pseudo-representative
                                        # out of it.
                                        a := ShallowCopy(reps[j][x]);
                                        a[3] := j;
                                    else
                                        a := reps[j][x];
                                    fi;
                                    a[4] := max1;
                                fi;
                                if  Length(y) <> 4  then
                                    if  y[2] <> k  then
                                        # we have to ensure that num(<y>) = k
                                        # when we construct a new
                                        # pseudo-representative out of it.
                                        b := ShallowCopy(y);
                                        b[2] := k;
                                    else
                                        b := y;
                                    fi;
                                    b[5] := max2;
                                else
                                    if  y[3] <> k  then
                                        # we have to ensure that num(<y>) = k
                                        # when we construct a new
                                        # pseudo-representative out of it.
                                        b := ShallowCopy(y);
                                        b[3] := k;
                                    else
                                        b := y;
                                    fi;
                                    b[4] := max2;
                                fi;
                                # now finally construct the
                                # pseudo-representative and add it to
                                # reps
                                z := [a, b, i, 0];
                                if  i >= avec[n]  then
                                    Add(reps[i], z);
                                else
                                    l := 3;
                                    while  l <= Length(pr[j][k])  and
                                      pr[j][k][l] < avec[n]  do
                                        Add(reps[ pr[j][k][l] ], z);
                                        l := l+2;
                                    od;
                                fi;
                            fi;
                        od;
                    od;
                fi;
                k := k+1;
            od;
        od;
    od;
    # now use the pseudo-representatives to get the desired polynomials
    for  i in [n..Length(pr)]  do
        for  j in [2..Length(reps[i])]  do
            # the first case: reps[i][j] is a "real" pseudo-representative
            if  Length(reps[i][j]) = 4  then
                if  reps[i][j][3] = i  then
                    GetPols(reps[i][j], pr, pols);
                fi;
            # the second case: reps[i][j] is a tree
            elif  reps[i][j][1] <> 0  then
                if  reps[i][j][2] = i  then
                    UnmarkTree(reps[i][j]);
                    hilf := MakeFormulaVector(reps[i][j], pr);
                    Dt_add(hilf, pols, pr);
                fi;
            # the third case: reps[i][j] is a deep thought monomial
            else
                Dt_add(reps[i][j], pols, pr);
            fi;
       od;
       Unbind(reps[i]);
   od;
   # finally convert the polynomials to the final state
   pols := Dt_Convert(pols);
   if  pols <> 0  then
      Dt_ConvertCoeffVecs(pols.evlistvec);
   fi;
   return pols;
end );


#############################################################################
##
#F  Calcreps2( <pr> ) . . . . . . . . . . compute the Deep-Thought-polynomials
##
##  'Calcreps2' returns the polynomials which have to be evaluated when
##  computing word*g_n^(y_n) for all <dtbound> <= n <= m where m is the
##  number of generators in the given presentation <pr>.
##
BindGlobal( "Calcreps2", function(pr, max, dtbound)
    local  i,reps,avec,max2, max1;

    reps := [];
    avec := Dt_Mkavec(pr);
    if  max >= Length(pr)  then
        max1 := Length(pr);
    else
        max1 := max;
    fi;
    for  i in [dtbound..Length(pr)]  do
        if  i >= max1  then
            max1 := Length(pr);
        fi;
        reps[i] := Calcrepsn(i, avec, pr, max1);
    od;
    max2 := 1;
    for  i in [1..Length(reps)]  do
        if  IsRecord(reps[i])  then
            max2 := i;
        fi;
    od;
    for  i in [1..max2]  do
        if  not IsRecord(reps[i])  then
            reps[i] := 1;
        fi;
    od;
    return reps;
end );

[ Dauer der Verarbeitung: 0.8 Sekunden  (vorverarbeitet)  ]