Quelle permutat.cc Sprache: C

/****************************************************************************
**
** This file is part of GAP, a system for computational discrete algebra.
**
** 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 contains the functions for permutations (small and large).
**
** Mathematically a permutation is a bijective mapping of a finite set onto
** itself. In \GAP\ this subset must always be of the form [ 1, 2, .., N ],
** where N is at most $2^32$.
**
** Internally a permutation is viewed as a mapping of [ 0, 1, .., N-1 ],
** because in C indexing of arrays is done with the origin 0 instead of 1.
** A permutation is represented by a bag of type 'T_PERM2' or 'T_PERM4' of
** the form
**
** (CONST_)ADDR_OBJ(p)
** |
** v
** +------+-------+-------+-------+-------+- - - -+-------+-------+
** | inv. | image | image | image | image | | image | image |
** | perm | of 0 | of 1 | of 2 | of 3 | | of N-2| of N-1|
** +------+-------+-------+-------+-------+- - - -+-------+-------+
** ^
** |
** (CONST_)ADDR_PERM2(p) resp. (CONST_)ADDR_PERM4(p)
**
** The first entry of the bag is either zero, or a reference to another
** permutation representing the inverse of . The remaining entries of the
** bag form an array describing the permutation. For bags of type 'T_PERM2',
** the entries are of type 'UInt2' (defined in 'common.h' as an 16 bit wide
** unsigned integer type), for type 'T_PERM4' the entries are of type
** 'UInt4' (defined as a 32bit wide unsigned integer type). The first of
** these entries is the image of 0, the second is the image of 1, and so on.
** Thus, the entry at C index is the image of , if we view the
** permutation as mapping of [ 0, 1, 2, .., N-1 ] as described above.
**
** Permutations are never shortened. For example, if the product of two
** permutations of degree 100 is the identity, it is nevertheless stored as
** array of length 100, in which the -th entry is of course simply .
** Testing whether a product has trailing fixpoints would be pretty costly,
** and permutations of equal degree can be handled by the functions faster.
*/

extern "C" {

#include "permutat.h"

#include "ariths.h"
#include "bool.h"
#include "error.h"
#include "gapstate.h"
#include "integer.h"
#include "io.h"
#include "listfunc.h"
#include "lists.h"
#include "modules.h"
#include "opers.h"
#include "plist.h"
#include "precord.h"
#include "range.h"
#include "records.h"
#include "saveload.h"
#include "sysfiles.h"
#include "trans.h"

} // extern "C"

#include "permutat_intern.hh"

/****************************************************************************
**
*V IdentityPerm . . . . . . . . . . . . . . . . . . . identity permutation
**
** 'IdentityPerm' is an identity permutation.
*/
Obj IdentityPerm;

static const UInt MAX_DEG_PERM4 = ((Int)1 << (sizeof(UInt) == 8 ? 32 : 28)) - 1;

static ModuleStateOffset PermutatStateOffset = -1;

typedef struct {

/****************************************************************************
**
*V TmpPerm . . . . . . . handle of the buffer bag of the permutation package
**
** 'TmpPerm' is the handle of a bag of type 'T_PERM4', which is created at
** initialization time of this package. Functions in this package can use
** this bag for whatever purpose they want. They have to make sure of
** course that it is large enough.
** The buffer is *not* guaranteed to have any particular value, routines
** that require a zero-initialization need to do this at the start.
** This buffer is only constructed once it is needed, to avoid startup
** costs (particularly when starting new threads).
** Use the UseTmpPerm(<size>) utility function to ensure it is constructed!
*/
Obj TmpPerm;

} PermutatModuleState;

#define TmpPerm MODULE_STATE(Permutat).TmpPerm

static UInt1 * UseTmpPerm(UInt size)
{
 if (TmpPerm == (Obj)0)
 TmpPerm = NewBag(T_PERM4, size);
 else if (SIZE_BAG(TmpPerm) < size)
 ResizeBag(TmpPerm, size);
 return (UInt1 *)(ADDR_OBJ(TmpPerm) + 1);
}

template <typename T>
static inline T * ADDR_TMP_PERM()
{
 // no GAP_ASSERT here on purpose
 return (T *)(ADDR_OBJ(TmpPerm) + 1);
}

/****************************************************************************
**
*F TypePerm( <perm> ) . . . . . . . . . . . . . . . . type of a permutation
**
** 'TypePerm' returns the type of permutations.
**
** 'TypePerm' is the function in 'TypeObjFuncs' for permutations.
*/
static Obj TYPE_PERM2;

static Obj TYPE_PERM4;

static Obj TypePerm2(Obj perm)
{
 return TYPE_PERM2;
}

static Obj TypePerm4(Obj perm)
{
 return TYPE_PERM4;
}

// Forward declarations
template <typename T>
static inline UInt LargestMovedPointPerm_(Obj perm);

template <typename T>
static Obj InvPerm(Obj perm);

/****************************************************************************
**
*F PrintPerm( <perm> ) . . . . . . . . . . . . . . . . . print a permutation
**
** 'PrintPerm' prints the permutation <perm> in the usual cycle notation. It
** uses the degree to print all points with same width, which looks nicer.
** Linebreaks are preferred most after cycles and next most after commas.
**
** It remembers which points have already been printed, to avoid O(n^2)
** behaviour.
*/
template <typename T>
static void PrintPerm(Obj perm)
{
 UInt degPerm; // degree of the permutation
 const T * ptPerm; // pointer to the permutation
 UInt p, q; // loop variables
 BOOL isId; // permutation is the identity?
 const char * fmt1; // common formats to print points
 const char * fmt2; // common formats to print points

 // set up the formats used, so all points are printed with equal width
 // Use LargestMovedPoint so printing is consistent
 degPerm = LargestMovedPointPerm_<T>(perm);
 if ( degPerm < 10 ) { fmt1 = "%>(%>%1d%<"; fmt2 = ",%>%1d%<"; }
 else if ( degPerm < 100 ) { fmt1 = "%>(%>%2d%<"; fmt2 = ",%>%2d%<"; }
 else if ( degPerm < 1000 ) { fmt1 = "%>(%>%3d%<"; fmt2 = ",%>%3d%<"; }
 else if ( degPerm < 10000 ) { fmt1 = "%>(%>%4d%<"; fmt2 = ",%>%4d%<"; }
 else { fmt1 = "%>(%>%5d%<"; fmt2 = ",%>%5d%<"; }

 UseTmpPerm(SIZE_OBJ(perm));
 T * ptSeen = ADDR_TMP_PERM<T>();

 // clear the buffer bag
 memset(ptSeen, 0, DEG_PERM<T>(perm) * sizeof(T));

 // run through all points
 isId = TRUE;
 ptPerm = CONST_ADDR_PERM<T>(perm);
 for ( p = 0; p < degPerm; p++ ) {
 // if the smallest is the one we started with lets print the cycle
 if (!ptSeen[p] && ptPerm[p] != p) {
 ptSeen[p] = 1;
 isId = FALSE;
 Pr(fmt1,(Int)(p+1), 0);
 // restore pointer, in case Pr caused a garbage collection
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptSeen = ADDR_TMP_PERM<T>();
 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
 ptSeen[q] = 1;
 Pr(fmt2,(Int)(q+1), 0);
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptSeen = ADDR_TMP_PERM<T>();
 }
 Pr("%<)", 0, 0);
 // restore pointer, in case Pr caused a garbage collection
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptSeen = ADDR_TMP_PERM<T>();
 }
 }

 // special case for the identity
 if ( isId ) Pr("()", 0, 0);
}

/****************************************************************************
**
*F EqPerm( <opL>, <opR> ) . . . . . . . test if two permutations are equal
**
** 'EqPerm' returns 'true' if the two permutations <opL> and <opR> are equal
** and 'false' otherwise.
**
** Two permutations may be equal, even if the two sequences do not have the
** same length, if the larger permutation fixes the exceeding points.
**
*/
template <typename TL, typename TR>
static Int EqPerm(Obj opL, Obj opR)
{
 UInt degL; // degree of the left operand
 const TL * ptL; // pointer to the left operand
 UInt degR; // degree of the right operand
 const TR * ptR; // pointer to the right operand
 UInt p; // loop variable

 // get the degrees
 degL = DEG_PERM<TL>(opL);
 degR = DEG_PERM<TR>(opR);

 // set up the pointers
 ptL = CONST_ADDR_PERM<TL>(opL);
 ptR = CONST_ADDR_PERM<TR>(opR);

 // search for a difference and return False if you find one
 if ( degL <= degR ) {
 for ( p = 0; p < degL; p++ )
 if ( *(ptL++) != *(ptR++) )
 return 0;
 for ( p = degL; p < degR; p++ )
 if ( p != *(ptR++) )
 return 0;
 }
 else {
 for ( p = 0; p < degR; p++ )
 if ( *(ptL++) != *(ptR++) )
 return 0;
 for ( p = degR; p < degL; p++ )
 if ( *(ptL++) != p )
 return 0;
 }

 // otherwise they must be equal
 return 1;
}

/****************************************************************************
**
*F LtPerm( <opL>, <opR> ) . test if one permutation is smaller than another
**
** 'LtPerm' returns 'true' if the permutation <opL> is strictly less than
** the permutation <opR>. Permutations are ordered lexicographically with
** respect to the images of 1,2,.., etc.
*/
template <typename TL, typename TR>
static Int LtPerm(Obj opL, Obj opR)
{
 UInt degL; // degree of the left operand
 const TL * ptL; // pointer to the left operand
 UInt degR; // degree of the right operand
 const TR * ptR; // pointer to the right operand
 UInt p; // loop variable

 // get the degrees of the permutations
 degL = DEG_PERM<TL>(opL);
 degR = DEG_PERM<TR>(opR);

 // set up the pointers
 ptL = CONST_ADDR_PERM<TL>(opL);
 ptR = CONST_ADDR_PERM<TR>(opR);

 // search for a difference and return if you find one
 if ( degL <= degR ) {

 for (p = 0; p < degL; p++, ptL++, ptR++)
 if (*ptL != *ptR) {
 return *ptL < *ptR;
 }
 for (p = degL; p < degR; p++, ptR++)
 if (p != *ptR) {
 return p < *ptR;
 }
 }
 else {
 for (p = 0; p < degR; p++, ptL++, ptR++)
 if (*ptL != *ptR) {
 return *ptL < *ptR;
 }
 for (p = degR; p < degL; p++, ptL++)
 if (*ptL != p) {
 return *ptL < p;
 }
 }

 // otherwise they must be equal
 return 0;
}

/****************************************************************************
**
*F ProdPerm( <opL>, <opR> ) . . . . . . . . . . . . product of permutations
**
** 'ProdPerm' returns the product of the two permutations <opL> and <opR>.
**
** This is a little bit tuned but should be sufficiently easy to understand.
*/
template <typename TL, typename TR>
static Obj ProdPerm(Obj opL, Obj opR)
{
 typedef typename ResultType<TL,TR>::type Res;

 Obj prd; // handle of the product (result)
 UInt degP; // degree of the product
 Res * ptP; // pointer to the product
 UInt degL; // degree of the left operand
 const TL * ptL; // pointer to the left operand
 UInt degR; // degree of the right operand
 const TR * ptR; // pointer to the right operand
 UInt p; // loop variable

 degL = DEG_PERM<TL>(opL);
 degR = DEG_PERM<TR>(opR);

 // handle trivial cases
 if (degL == 0) {
 return opR;
 }
 if (degR == 0) {
 return opL;
 }

 // compute the size of the result and allocate a bag
 degP = degL < degR ? degR : degL;
 prd = NEW_PERM<Res>( degP );

 // set up the pointers
 ptL = CONST_ADDR_PERM<TL>(opL);
 ptR = CONST_ADDR_PERM<TR>(opR);
 ptP = ADDR_PERM<Res>(prd);

 // if the left (inner) permutation has smaller degree, it is very easy
 if ( degL <= degR ) {
 for ( p = 0; p < degL; p++ )
 *(ptP++) = ptR[ *(ptL++) ];
 for ( p = degL; p < degR; p++ )
 *(ptP++) = ptR[ p ];
 }

 // otherwise we have to use the macro 'IMAGE'
 else {
 for ( p = 0; p < degL; p++ )
 *(ptP++) = IMAGE( ptL[ p ], ptR, degR );
 }

 return prd;
}

/****************************************************************************
**
*F QuoPerm( <opL>, <opR> ) . . . . . . . . . . . . quotient of permutations
**
** 'QuoPerm' returns the quotient of the permutations <opL> and <opR>, i.e.,
** the product '<opL>\*<opR>\^-1'.
**
** Unfortunately this cannot be done in <degree> steps, we need 2 * <degree>
** steps.
*/
static Obj QuoPerm(Obj opL, Obj opR)
{
 return PROD(opL, INV(opR));
}

/****************************************************************************
**
*F LQuoPerm( <opL>, <opR> ) . . . . . . . . . left quotient of permutations
**
** 'LQuoPerm' returns the left quotient of the two permutations <opL> and
** <opR>, i.e., the value of '<opL>\^-1*<opR>', which sometimes comes handy.
**
** This can be done as fast as a single multiplication or inversion.
*/
template <typename TL, typename TR>
static Obj LQuoPerm(Obj opL, Obj opR)
{
 typedef typename ResultType<TL,TR>::type Res;

 Obj mod; // handle of the quotient (result)
 UInt degM; // degree of the quotient
 Res * ptM; // pointer to the quotient
 UInt degL; // degree of the left operand
 const TL * ptL; // pointer to the left operand
 UInt degR; // degree of the right operand
 const TR * ptR; // pointer to the right operand
 UInt p; // loop variable

 degL = DEG_PERM<TL>(opL);
 degR = DEG_PERM<TR>(opR);

 // handle trivial cases
 if (degL == 0) {
 return opR;
 }
 if (degR == 0) {
 return InvPerm<TL>(opL);
 }

 // compute the size of the result and allocate a bag
 degM = degL < degR ? degR : degL;
 mod = NEW_PERM<Res>( degM );

 // set up the pointers
 ptL = CONST_ADDR_PERM<TL>(opL);
 ptR = CONST_ADDR_PERM<TR>(opR);
 ptM = ADDR_PERM<Res>(mod);

 // it is one thing if the left (inner) permutation is smaller
 if ( degL <= degR ) {
 for ( p = 0; p < degL; p++ )
 ptM[ *(ptL++) ] = *(ptR++);
 for ( p = degL; p < degR; p++ )
 ptM[ p ] = *(ptR++);
 }

 // and another if the right (outer) permutation is smaller
 else {
 for ( p = 0; p < degR; p++ )
 ptM[ *(ptL++) ] = *(ptR++);
 for ( p = degR; p < degL; p++ )
 ptM[ *(ptL++) ] = p;
 }

 return mod;
}

/****************************************************************************
**
*F InvPerm( <perm> ) . . . . . . . . . . . . . . . inverse of a permutation
*/
template <typename T>
static Obj InvPerm(Obj perm)
{
 Obj inv; // handle of the inverse (result)
 T * ptInv; // pointer to the inverse
 const T * ptPerm; // pointer to the permutation
 UInt deg; // degree of the permutation
 UInt p; // loop variables

 inv = STOREDINV_PERM(perm);
 if (inv != 0)
 return inv;

 deg = DEG_PERM<T>(perm);
 inv = NEW_PERM<T>(deg);

 // get pointer to the permutation and the inverse
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptInv = ADDR_PERM<T>(inv);

 // invert the permutation
 for ( p = 0; p < deg; p++ )
 ptInv[ *ptPerm++ ] = p;

 // store and return the inverse
 SET_STOREDINV_PERM(perm, inv);
 return inv;
}

/****************************************************************************
**
*F PowPermInt( <opL>, <opR> ) . . . . . . . integer power of a permutation
**
** 'PowPermInt' returns the <opR>-th power of the permutation <opL>. <opR>
** must be a small integer.
**
** This repeatedly applies the permutation <opR> to all points which seems
** to be faster than binary powering, and does not need temporary storage.
*/
template <typename T>
static Obj PowPermInt(Obj opL, Obj opR)
{
 Obj pow; // handle of the power (result)
 T * ptP; // pointer to the power
 const T * ptL; // pointer to the permutation
 UInt1 * ptKnown; // pointer to temporary bag
 UInt deg; // degree of the permutation
 Int exp, e; // exponent (right operand)
 UInt len; // length of cycle (result)
 UInt p, q, r; // loop variables

 // handle zeroth and first powers and stored inverses separately
 if ( opR == INTOBJ_INT(0))
 return IdentityPerm;
 if ( opR == INTOBJ_INT(1))
 return opL;
 if (opR == INTOBJ_INT(-1))
 return InvPerm<T>(opL);

 deg = DEG_PERM<T>(opL);

 // handle trivial permutations
 if (deg == 0) {
 return IdentityPerm;
 }

 // allocate a result bag
 pow = NEW_PERM<T>(deg);

 // compute the power by repeated mapping for small positive exponents
 if ( IS_INTOBJ(opR)
 && 2 <= INT_INTOBJ(opR) && INT_INTOBJ(opR) < 8 ) {

 // get pointer to the permutation and the power
 exp = INT_INTOBJ(opR);
 ptL = CONST_ADDR_PERM<T>(opL);
 ptP = ADDR_PERM<T>(pow);

 // loop over the points of the permutation
 for ( p = 0; p < deg; p++ ) {
 q = p;
 for ( e = 0; e < exp; e++ )
 q = ptL[q];
 ptP[p] = q;
 }

 }

 // compute the power by raising the cycles individually for large exps
 else if ( IS_INTOBJ(opR) && 8 <= INT_INTOBJ(opR) ) {

 // make sure that the buffer bag is large enough
 UseTmpPerm(SIZE_OBJ(opL));
 ptKnown = ADDR_TMP_PERM<UInt1>();

 // clear the buffer bag
 memset(ptKnown, 0, DEG_PERM<T>(opL));

 // get pointer to the permutation and the power
 exp = INT_INTOBJ(opR);
 ptL = CONST_ADDR_PERM<T>(opL);
 ptP = ADDR_PERM<T>(pow);

 // loop over all cycles
 for ( p = 0; p < deg; p++ ) {

 // if we haven't looked at this cycle so far
 if ( ptKnown[p] == 0 ) {

 // find the length of this cycle
 len = 1;
 for ( q = ptL[p]; q != p; q = ptL[q] ) {
 len++; ptKnown[q] = 1;
 }

 // raise this cycle to the power <exp> mod <len>
 r = p;
 for ( e = 0; e < exp % len; e++ )
 r = ptL[r];
 ptP[p] = r;
 r = ptL[r];
 for ( q = ptL[p]; q != p; q = ptL[q] ) {
 ptP[q] = r;
 r = ptL[r];
 }

 }

 }

 }

 // compute the power by raising the cycles individually for large exps
 else if ( TNUM_OBJ(opR) == T_INTPOS ) {

 // make sure that the buffer bag is large enough
 UseTmpPerm(SIZE_OBJ(opL));
 ptKnown = ADDR_TMP_PERM<UInt1>();

 // clear the buffer bag
 memset(ptKnown, 0, DEG_PERM<T>(opL));

 // get pointer to the permutation and the power
 ptL = CONST_ADDR_PERM<T>(opL);
 ptP = ADDR_PERM<T>(pow);

 // loop over all cycles
 for ( p = 0; p < deg; p++ ) {

 // if we haven't looked at this cycle so far
 if ( ptKnown[p] == 0 ) {

 // find the length of this cycle
 len = 1;
 for ( q = ptL[p]; q != p; q = ptL[q] ) {
 len++; ptKnown[q] = 1;
 }

 // raise this cycle to the power <exp> mod <len>
 r = p;
 exp = INT_INTOBJ( ModInt( opR, INTOBJ_INT(len) ) );
 for ( e = 0; e < exp; e++ )
 r = ptL[r];
 ptP[p] = r;
 r = ptL[r];
 for ( q = ptL[p]; q != p; q = ptL[q] ) {
 ptP[q] = r;
 r = ptL[r];
 }

 }

 }

 }

 // compute the power by repeated mapping for small negative exponents
 else if ( IS_INTOBJ(opR)
 && -8 < INT_INTOBJ(opR) && INT_INTOBJ(opR) < 0 ) {

 // get pointer to the permutation and the power
 exp = -INT_INTOBJ(opR);
 ptL = CONST_ADDR_PERM<T>(opL);
 ptP = ADDR_PERM<T>(pow);

 // loop over the points
 for ( p = 0; p < deg; p++ ) {
 q = p;
 for ( e = 0; e < exp; e++ )
 q = ptL[q];
 ptP[q] = p;
 }

 }

 // compute the power by raising the cycles individually for large exps
 else if ( IS_INTOBJ(opR) && INT_INTOBJ(opR) <= -8 ) {

 // make sure that the buffer bag is large enough
 UseTmpPerm(SIZE_OBJ(opL));
 ptKnown = ADDR_TMP_PERM<UInt1>();

 // clear the buffer bag
 memset(ptKnown, 0, DEG_PERM<T>(opL));

 // get pointer to the permutation and the power
 exp = -INT_INTOBJ(opR);
 ptL = CONST_ADDR_PERM<T>(opL);
 ptP = ADDR_PERM<T>(pow);

 // loop over all cycles
 for ( p = 0; p < deg; p++ ) {

 // if we haven't looked at this cycle so far
 if ( ptKnown[p] == 0 ) {

 // find the length of this cycle
 len = 1;
 for ( q = ptL[p]; q != p; q = ptL[q] ) {
 len++; ptKnown[q] = 1;
 }

 // raise this cycle to the power <exp> mod <len>
 r = p;
 for ( e = 0; e < exp % len; e++ )
 r = ptL[r];
 ptP[r] = p;
 r = ptL[r];
 for ( q = ptL[p]; q != p; q = ptL[q] ) {
 ptP[r] = q;
 r = ptL[r];
 }

 }

 }

 }

 // compute the power by raising the cycles individually for large exps
 else if ( TNUM_OBJ(opR) == T_INTNEG ) {
 // do negation first as it can cause a garbage collection
 opR = AInvInt(opR);

 // make sure that the buffer bag is large enough
 UseTmpPerm(SIZE_OBJ(opL));
 ptKnown = ADDR_TMP_PERM<UInt1>();

 // clear the buffer bag
 memset(ptKnown, 0, DEG_PERM<T>(opL));

 // get pointer to the permutation and the power
 ptL = CONST_ADDR_PERM<T>(opL);
 ptP = ADDR_PERM<T>(pow);

 // loop over all cycles
 for ( p = 0; p < deg; p++ ) {

 // if we haven't looked at this cycle so far
 if ( ptKnown[p] == 0 ) {

 // find the length of this cycle
 len = 1;
 for ( q = ptL[p]; q != p; q = ptL[q] ) {
 len++; ptKnown[q] = 1;
 }

 // raise this cycle to the power <exp> mod <len>
 r = p;
 exp = INT_INTOBJ( ModInt( opR, INTOBJ_INT(len) ) );
 for ( e = 0; e < exp % len; e++ )
 r = ptL[r];
 ptP[r] = p;
 r = ptL[r];
 for ( q = ptL[p]; q != p; q = ptL[q] ) {
 ptP[r] = q;
 r = ptL[r];
 }

 }

 }

 }

 return pow;
}

/****************************************************************************
**
*F PowIntPerm( <opL>, <opR> ) . . . image of an integer under a permutation
**
** 'PowIntPerm' returns the image of the positive integer <opL> under the
** permutation <opR>. If <opL> is larger than the degree of <opR> it is a
** fixpoint of the permutation and thus simply returned.
*/
template <typename T>
static Obj PowIntPerm(Obj opL, Obj opR)
{
 Int img; // image (result)

 GAP_ASSERT(TNUM_OBJ(opL) == T_INTPOS || TNUM_OBJ(opL) == T_INT);

 // large positive integers (> 2^28-1) are fixed by any permutation
 if ( TNUM_OBJ(opL) == T_INTPOS )
 return opL;

 img = INT_INTOBJ(opL);
 RequireArgumentConditionEx("PowIntPerm", opL, "", img > 0,
 "must be a positive integer");

 // compute the image
 if ( img <= DEG_PERM<T>(opR) ) {
 img = (CONST_ADDR_PERM<T>(opR))[img-1] + 1;
 }

 // return it
 return INTOBJ_INT(img);
}

/****************************************************************************
**
*F QuoIntPerm( <opL>, <opR> ) . preimage of an integer under a permutation
**
** 'QuoIntPerm' returns the preimage of the preimage integer <opL> under the
** permutation <opR>. If <opL> is larger than the degree of <opR> it is a
** fixpoint, and thus simply returned.
**
** There are basically two ways to find the preimage. One is to run through
** <opR> and look for <opL>. The index where it's found is the preimage.
** The other is to find the image of <opL> under <opR>, the image of that
** point and so on, until we come back to <opL>. The last point is the
** preimage of <opL>. This is faster because the cycles are usually short.
*/
static Obj PERM_INVERSE_THRESHOLD;

template <typename T>
static Obj QuoIntPerm(Obj opL, Obj opR)
{
 T pre; // preimage (result)
 Int img; // image (left operand)
 const T * ptR; // pointer to the permutation

 GAP_ASSERT(TNUM_OBJ(opL) == T_INTPOS || TNUM_OBJ(opL) == T_INT);

 // large positive integers (> 2^28-1) are fixed by any permutation
 if ( TNUM_OBJ(opL) == T_INTPOS )
 return opL;

 img = INT_INTOBJ(opL);
 RequireArgumentConditionEx("QuoIntPerm", opL, "", img > 0,
 "must be a positive integer");

 Obj inv = STOREDINV_PERM(opR);

 if (inv == 0 && PERM_INVERSE_THRESHOLD != 0 &&
 IS_INTOBJ(PERM_INVERSE_THRESHOLD) &&
 DEG_PERM<T>(opR) <= INT_INTOBJ(PERM_INVERSE_THRESHOLD))
 inv = InvPerm<T>(opR);

 if (inv != 0)
 return INTOBJ_INT(
 IMAGE(img - 1, CONST_ADDR_PERM<T>(inv), DEG_PERM<T>(inv)) + 1);

 // compute the preimage
 if ( img <= DEG_PERM<T>(opR) ) {
 pre = T(img - 1);
 ptR = CONST_ADDR_PERM<T>(opR);
 while (ptR[pre] != T(img - 1))
 pre = ptR[pre];
 // return it
 return INTOBJ_INT(pre + 1);
 }
 else
 return INTOBJ_INT(img);
}

/****************************************************************************
**
*F PowPerm( <opL>, <opR> ) . . . . . . . . . . . conjugation of permutations
**
** 'PowPerm' returns the conjugation of the two permutations <opL> and
** <opR>, that s defined as the following product '<opR>\^-1 \*\ <opL> \*\
** <opR>'.
*/
template <typename TL, typename TR>
static Obj PowPerm(Obj opL, Obj opR)
{
 typedef typename ResultType<TL,TR>::type Res;

 Obj cnj; // handle of the conjugation (res)
 UInt degC; // degree of the conjugation
 Res * ptC; // pointer to the conjugation
 UInt degL; // degree of the left operand
 const TL * ptL; // pointer to the left operand
 UInt degR; // degree of the right operand
 const TR * ptR; // pointer to the right operand
 UInt p; // loop variable

 degL = DEG_PERM<TL>(opL);
 degR = DEG_PERM<TR>(opR);

 // handle trivial cases
 if (degL == 0) {
 return IdentityPerm;
 }

 if (degR == 0) {
 return opL;
 }

 // compute the size of the result and allocate a bag
 degC = degL < degR ? degR : degL;
 cnj = NEW_PERM<Res>( degC );

 // set up the pointers
 ptL = CONST_ADDR_PERM<TL>(opL);
 ptR = CONST_ADDR_PERM<TR>(opR);
 ptC = ADDR_PERM<Res>(cnj);

 // it is faster if the both permutations have the same size
 if ( degL == degR ) {
 for ( p = 0; p < degC; p++ )
 ptC[ ptR[p] ] = ptR[ ptL[p] ];
 }

 // otherwise we have to use the macro 'IMAGE' three times
 else {
 for ( p = 0; p < degC; p++ )
 ptC[ IMAGE(p,ptR,degR) ] = IMAGE( IMAGE(p,ptL,degL), ptR, degR );
 }

 return cnj;
}

/****************************************************************************
**
*F CommPerm( <opL>, <opR> ) . . . . . . . . commutator of two permutations
**
** 'CommPerm' returns the commutator of the two permutations <opL> and
** <opR>, that is defined as '<hd>\^-1 \*\ <opR>\^-1 \*\ <opL> \*\ <opR>'.
*/
template <typename TL, typename TR>
static Obj CommPerm(Obj opL, Obj opR)
{
 typedef typename ResultType<TL,TR>::type Res;

 Obj com; // handle of the commutator (res)
 UInt degC; // degree of the commutator
 Res * ptC; // pointer to the commutator
 UInt degL; // degree of the left operand
 const TL * ptL; // pointer to the left operand
 UInt degR; // degree of the right operand
 const TR * ptR; // pointer to the right operand
 UInt p; // loop variable

 degL = DEG_PERM<TL>(opL);
 degR = DEG_PERM<TR>(opR);

 // handle trivial cases
 if (degL == 0 || degR == 0) {
 return IdentityPerm;
 }

 // compute the size of the result and allocate a bag
 degC = degL < degR ? degR : degL;
 com = NEW_PERM<Res>( degC );

 // set up the pointers
 ptL = CONST_ADDR_PERM<TL>(opL);
 ptR = CONST_ADDR_PERM<TR>(opR);
 ptC = ADDR_PERM<Res>(com);

 // it is faster if the both permutations have the same size
 if ( degL == degR ) {
 for ( p = 0; p < degC; p++ )
 ptC[ ptL[ ptR[ p ] ] ] = ptR[ ptL[ p ] ];
 }

 // otherwise we have to use the macro 'IMAGE' four times
 else {
 for ( p = 0; p < degC; p++ )
 ptC[ IMAGE( IMAGE(p,ptR,degR), ptL, degL ) ]
 = IMAGE( IMAGE(p,ptL,degL), ptR, degR );
 }

 return com;
}

/****************************************************************************
**
*F OnePerm( <perm> )
*/
static Obj OnePerm(Obj op)
{
 return IdentityPerm;
}

/****************************************************************************
**
*F FiltIS_PERM( <self>, <val> ) . . . . . . test if a value is a permutation
**
** 'FiltIS_PERM' implements the internal function 'IsPerm'.
**
** 'IsPerm( <val> )'
**
** 'IsPerm' returns 'true' if the value <val> is a permutation and 'false'
** otherwise.
*/
static Obj IsPermFilt;

static Obj FiltIS_PERM(Obj self, Obj val)
{
 // return 'true' if <val> is a permutation and 'false' otherwise
 if ( TNUM_OBJ(val) == T_PERM2 || TNUM_OBJ(val) == T_PERM4 ) {
 return True;
 }
 else if ( TNUM_OBJ(val) < FIRST_EXTERNAL_TNUM ) {
 return False;
 }
 else {
 return DoFilter( self, val );
 }
}

/****************************************************************************
**
*F FuncPermList( <self>, <list> ) . . . . . convert a list to a permutation
**
** 'FuncPermList' implements the internal function 'PermList'
**
** 'PermList( <list> )'
**
** Converts the list <list> into a permutation, which is then returned.
**
** 'FuncPermList' simply copies the list pointwise into a permutation bag.
** It also does some checks to make sure that the list is a permutation.
*/
template <typename T>
static inline Obj PermList(Obj list)
{
 Obj perm; // handle of the permutation
 T * ptPerm; // pointer to the permutation
 UInt degPerm; // degree of the permutation
 const Obj * ptList; // pointer to the list
 T * ptTmp; // pointer to the buffer bag
 Int i, k; // loop variables

 GAP_ASSERT(IS_PLIST(list));
 degPerm = LEN_PLIST( list );

 /* make sure that the global buffer bag is large enough for checkin*/
 UseTmpPerm(SIZEBAG_PERM<T>(degPerm));

 // allocate the bag for the permutation and get pointer
 perm = NEW_PERM<T>(degPerm);
 ptPerm = ADDR_PERM<T>(perm);
 ptList = CONST_ADDR_OBJ(list);
 ptTmp = ADDR_TMP_PERM<T>();

 // make the buffer bag clean
 memset(ptTmp, 0, degPerm * sizeof(T));

 // run through all entries of the list
 for ( i = 1; i <= degPerm; i++ ) {

 // get the th entry of the list
 if ( !IS_INTOBJ(ptList[i]) ) {
 return Fail;
 }
 k = INT_INTOBJ(ptList[i]);
 if ( k <= 0 || degPerm < k ) {
 return Fail;
 }

 // make sure we haven't seen this entry yet
 if ( ptTmp[k-1] != 0 ) {
 return Fail;
 }
 ptTmp[k-1] = 1;

 // and finally copy it into the permutation
 ptPerm[i-1] = k-1;
 }

 // return the permutation
 return perm;
}

static Obj FuncPermList(Obj self, Obj list)
{
 RequireSmallList(SELF_NAME, list);

 UInt len = LEN_LIST( list );
 if (len == 0)
 return IdentityPerm;
 if (!IS_PLIST(list)) {
 if (!IS_POSS_LIST(list))
 return Fail;
 if (IS_RANGE(list)) {
 if (GET_LOW_RANGE(list) == 1 && GET_INC_RANGE(list) == 1)
 return IdentityPerm;
 }
 list = PLAIN_LIST_COPY(list);
 }

 if ( len <= 65536 ) {
 return PermList<UInt2>(list);
 }
 else if (len <= MAX_DEG_PERM4) {
 return PermList<UInt4>(list);
 }
 else {
 ErrorMayQuit("PermList: list length %d exceeds maximum permutation degree",
 len, 0);
 }
}

template <typename T>
static inline Obj ListPerm_(Obj perm, Int len)
{
 Obj res; // handle of the image, result
 Obj * ptRes; // pointer to the result
 const T * ptPrm; // pointer to the permutation
 UInt deg; // degree of the permutation
 UInt i; // loop variable

 if (len <= 0)
 return NewEmptyPlist();

 // copy the list into a mutable plist, which we will then modify in place
 res = NEW_PLIST(T_PLIST_CYC, len);
 SET_LEN_PLIST(res, len);

 // get the pointer
 ptRes = ADDR_OBJ(res) + 1;
 ptPrm = CONST_ADDR_PERM<T>(perm);

 // loop over the entries of the permutation
 deg = DEG_PERM<T>(perm);
 if (deg > len)
 deg = len;
 for (i = 1; i <= deg; i++, ptRes++) {
 *ptRes = INTOBJ_INT(ptPrm[i - 1] + 1);
 }
 // add extras if requested
 for (; i <= len; i++, ptRes++) {
 *ptRes = INTOBJ_INT(i);
 }

 return res;
}

static Obj ListPermOper;

static Obj FuncListPerm1(Obj self, Obj perm)
{
 RequirePermutation(SELF_NAME, perm);
 Int nn = LargestMovedPointPerm(perm);
 if (TNUM_OBJ(perm) == T_PERM2)
 return ListPerm_<UInt2>(perm, nn);
 else
 return ListPerm_<UInt4>(perm, nn);
}

static Obj FuncListPerm2(Obj self, Obj perm, Obj n)
{
 RequirePermutation(SELF_NAME, perm);
 Int nn = GetSmallInt(SELF_NAME, n);
 if (TNUM_OBJ(perm) == T_PERM2)
 return ListPerm_<UInt2>(perm, nn);
 else
 return ListPerm_<UInt4>(perm, nn);
}

/****************************************************************************
**
*F LargestMovedPointPerm( <perm> ) largest point moved by perm
**
** 'LargestMovedPointPerm' returns the largest positive integer that is
** moved by the permutation <perm>.
**
** This is easy, except that permutations may contain trailing fixpoints.
*/
template <typename T>
static inline UInt LargestMovedPointPerm_(Obj perm)
{
 UInt sup;
 const T * ptPerm;

 ptPerm = CONST_ADDR_PERM<T>(perm);
 sup = DEG_PERM<T>(perm);
 while (sup-- > 0) {
 if (ptPerm[sup] != sup)
 return sup + 1;
 }
 return 0;
}

UInt LargestMovedPointPerm(Obj perm)
{
 GAP_ASSERT(TNUM_OBJ(perm) == T_PERM2 || TNUM_OBJ(perm) == T_PERM4);

 if (TNUM_OBJ(perm) == T_PERM2)
 return LargestMovedPointPerm_<UInt2>(perm);
 else
 return LargestMovedPointPerm_<UInt4>(perm);
}

/****************************************************************************
**
*F SmallestMovedPointPerm( <perm> ) smallest point moved by perm
**
** 'SmallestMovedPointPerm' implements 'SmallestMovedPoint' for internal
** permutations.
*/

// Import 'Infinity', as a return value for the identity permutation
static Obj Infinity;

template <typename T>
static inline Obj SmallestMovedPointPerm_(Obj perm)
{
 const T * ptPerm = CONST_ADDR_PERM<T>(perm);
 UInt deg = DEG_PERM<T>(perm);

 for (UInt i = 0; i < deg; i++) {
 if (ptPerm[i] != i)
 return INTOBJ_INT(i + 1);
 }
 return Infinity;
}

static Obj SmallestMovedPointPerm(Obj perm)
{
 GAP_ASSERT(TNUM_OBJ(perm) == T_PERM2 || TNUM_OBJ(perm) == T_PERM4);

 if (TNUM_OBJ(perm) == T_PERM2)
 return SmallestMovedPointPerm_<UInt2>(perm);
 else
 return SmallestMovedPointPerm_<UInt4>(perm);
}

/****************************************************************************
**
*F FuncLARGEST_MOVED_POINT_PERM( <self>, <perm> ) largest point moved by perm
**
** GAP-level wrapper for 'LargestMovedPointPerm'.
*/
static Obj FuncLARGEST_MOVED_POINT_PERM(Obj self, Obj perm)
{
 RequirePermutation(SELF_NAME, perm);

 return INTOBJ_INT(LargestMovedPointPerm(perm));
}

/****************************************************************************
**
*F FuncSMALLEST_MOVED_POINT_PERM( <self>, <perm> )
**
** GAP-level wrapper for 'SmallestMovedPointPerm'.
*/
static Obj FuncSMALLEST_MOVED_POINT_PERM(Obj self, Obj perm)
{
 RequirePermutation(SELF_NAME, perm);

 return SmallestMovedPointPerm(perm);
}

/****************************************************************************
**
*F FuncCYCLE_LENGTH_PERM_INT( <self>, <perm>, <point> ) . . . . . . . . . .
*F . . . . . . . . . . . . . . . . . . length of a cycle under a permutation
**
** 'FuncCycleLengthInt' implements the internal function
** 'CycleLengthPermInt'
**
** 'CycleLengthPermInt( <perm>, <point> )'
**
** 'CycleLengthPermInt' returns the length of the cycle of <point>, which
** must be a positive integer, under the permutation <perm>.
**
** Note that the order of the arguments to this function has been reversed.
*/
template <typename T>
static inline Obj CYCLE_LENGTH_PERM_INT(Obj perm, UInt pnt)
{
 const T * ptPerm; // pointer to the permutation
 UInt deg; // degree of the permutation
 UInt len; // length of cycle (result)
 UInt p; // loop variable

 // get pointer to the permutation, the degree, and the point
 ptPerm = CONST_ADDR_PERM<T>(perm);
 deg = DEG_PERM<T>(perm);

 // now compute the length by looping over the cycle
 len = 1;
 if ( pnt < deg ) {
 for ( p = ptPerm[pnt]; p != pnt; p = ptPerm[p] )
 len++;
 }

 return INTOBJ_INT(len);
}

static Obj FuncCYCLE_LENGTH_PERM_INT(Obj self, Obj perm, Obj point)
{
 UInt pnt; // value of the point

 RequirePermutation(SELF_NAME, perm);
 pnt = GetPositiveSmallInt("CycleLengthPermInt", point) - 1;

 if ( TNUM_OBJ(perm) == T_PERM2 ) {
 return CYCLE_LENGTH_PERM_INT<UInt2>(perm, pnt);
 }
 else {
 return CYCLE_LENGTH_PERM_INT<UInt4>(perm, pnt);
 }
}

/****************************************************************************
**
*F FuncCYCLE_PERM_INT( <self>, <perm>, <point> ) . . cycle of a permutation
*
** 'FuncCYCLE_PERM_INT' implements the internal function 'CyclePermInt'.
**
** 'CyclePermInt( <perm>, <point> )'
**
** 'CyclePermInt' returns the cycle of <point>, which must be a positive
** integer, under the permutation <perm> as a list.
*/
template <typename T>
static inline Obj CYCLE_PERM_INT(Obj perm, UInt pnt)
{
 Obj list; // handle of the list (result)
 Obj * ptList; // pointer to the list
 const T * ptPerm; // pointer to the permutation
 UInt deg; // degree of the permutation
 UInt len; // length of the cycle
 UInt p; // loop variable

 // get pointer to the permutation, the degree, and the point
 ptPerm = CONST_ADDR_PERM<T>(perm);
 deg = DEG_PERM<T>(perm);

 // now compute the length by looping over the cycle
 len = 1;
 if ( pnt < deg ) {
 for ( p = ptPerm[pnt]; p != pnt; p = ptPerm[p] )
 len++;
 }

 // allocate the list
 list = NEW_PLIST( T_PLIST, len );
 SET_LEN_PLIST( list, len );
 ptList = ADDR_OBJ(list);
 ptPerm = CONST_ADDR_PERM<T>(perm);

 // copy the points into the list
 len = 1;
 ptList[len++] = INTOBJ_INT( pnt+1 );
 if ( pnt < deg ) {
 for ( p = ptPerm[pnt]; p != pnt; p = ptPerm[p] )
 ptList[len++] = INTOBJ_INT( p+1 );
 }

 return list;
}

static Obj FuncCYCLE_PERM_INT(Obj self, Obj perm, Obj point)
{
 UInt pnt; // value of the point

 RequirePermutation(SELF_NAME, perm);
 pnt = GetPositiveSmallInt("CyclePermInt", point) - 1;

 if ( TNUM_OBJ(perm) == T_PERM2 ) {
 return CYCLE_PERM_INT<UInt2>(perm, pnt);
 }
 else {
 return CYCLE_PERM_INT<UInt4>(perm, pnt);
 }
}

/****************************************************************************
**
*F FuncCYCLE_STRUCT_PERM( <self>, <perm> ) . . . . . cycle structure of perm
*
** 'FuncCYCLE_STRUCT_PERM' implements the internal function
** `CycleStructPerm'.
**
** `CycleStructPerm( <perm> )'
**
** 'CycleStructPerm' returns a list of the form as described under
** `CycleStructure'.
*/
template <typename T>
static inline Obj CYCLE_STRUCT_PERM(Obj perm)
{
 Obj list; // handle of the list (result)
 Obj * ptList; // pointer to the list
 const T * ptPerm; // pointer to the permutation
 T * scratch;
 T * offset;
 UInt deg; // degree of the permutation
 UInt pnt; // value of the point
 UInt len; // length of the cycle
 UInt p; // loop variable
 UInt max; // maximal cycle length
 UInt cnt;
 UInt ende;
 UInt bytes;
 UInt1 * clr;

 // make sure that the buffer bag is large enough
 UseTmpPerm(SIZE_OBJ(perm) + 8);

 // find the largest moved point
 deg = LargestMovedPointPerm_<T>(perm);
 if (deg == 0) {
 // special treatment of identity
 return NewEmptyPlist();
 }

 scratch = ADDR_TMP_PERM<T>();

 // the first deg bytes of TmpPerm hold a bit list of points done
 // so far. The remaining bytes will form the lengths of nontrivial
 // cycles (as numbers of type T). As every nontrivial cycle requires
 // at least 2 points, this is guaranteed to fit.
 bytes = ((deg / sizeof(T)) + 1) * sizeof(T); // ensure alignment
 offset = (T *)((UInt)scratch + (bytes));

 // clear out the bits
 memset(scratch, 0, bytes);

 cnt = 0;
 clr = (UInt1 *)scratch;
 max = 0;
 ptPerm = CONST_ADDR_PERM<T>(perm);
 for (pnt = 0; pnt < deg; pnt++) {
 if (clr[pnt] == 0) {
 len = 0;
 clr[pnt] = 1;
 for (p = ptPerm[pnt]; p != pnt; p = ptPerm[p]) {
 clr[p] = 1;
 len++;
 }

 if (len > 0) {
 offset[cnt] = (T)len;
 cnt++;
 if (len > max) {
 max = len;
 }
 }
 }
 }

 ende = cnt;

 // create the list
 list = NEW_PLIST(T_PLIST, max);
 SET_LEN_PLIST(list, max);
 ptList = ADDR_OBJ(list);

 // Recalculate after possible GC
 scratch = ADDR_TMP_PERM<T>();
 offset = (T *)((UInt)scratch + (bytes));

 for (cnt = 0; cnt < ende; cnt++) {
 pnt = (UInt)offset[cnt];
 if (ptList[pnt] == 0) {
 ptList[pnt] = INTOBJ_INT(1);
 }
 else {
 ptList[pnt] = INTOBJ_INT(INT_INTOBJ(ptList[pnt]) + 1);
 }
 }

 return list;
}

static Obj FuncCYCLE_STRUCT_PERM(Obj self, Obj perm)
{
 RequirePermutation(SELF_NAME, perm);

 if (TNUM_OBJ(perm) == T_PERM2) {
 return CYCLE_STRUCT_PERM<UInt2>(perm);
 }
 else {
 return CYCLE_STRUCT_PERM<UInt4>(perm);
 }
}

/****************************************************************************
**
*F FuncORDER_PERM( <self>, <perm> ) . . . . . . . . . order of a permutation
**
** 'FuncORDER_PERM' implements the internal function 'OrderPerm'.
**
** 'OrderPerm( <perm> )'
**
** 'OrderPerm' returns the order of the permutation <perm>, i.e., the
** smallest positive integer <n> such that '<perm>\^<n>' is the identity.
**
** Since the largest element in S(65536) has order greater than 10^382 this
** computation may easily overflow. So we have to use arbitrary precision.
*/
template <typename T>
static inline Obj ORDER_PERM(Obj perm)
{
 const T * ptPerm; // pointer to the permutation
 Obj ord; // order (result), may be huge
 T * ptKnown; // pointer to temporary bag
 UInt len; // length of one cycle
 UInt p, q; // loop variables

 // make sure that the buffer bag is large enough
 UseTmpPerm(SIZE_OBJ(perm));

 // get the pointer to the bags
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptKnown = ADDR_TMP_PERM<T>();

 // clear the buffer bag
 for ( p = 0; p < DEG_PERM<T>(perm); p++ )
 ptKnown[p] = 0;

 // start with order 1
 ord = INTOBJ_INT(1);

 // loop over all cycles
 for ( p = 0; p < DEG_PERM<T>(perm); p++ ) {

 // if we haven't looked at this cycle so far
 if ( ptKnown[p] == 0 && ptPerm[p] != p ) {

 // find the length of this cycle
 len = 1;
 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
 len++; ptKnown[q] = 1;
 }

 ord = LcmInt( ord, INTOBJ_INT( len ) );

 // update bag pointers, in case a garbage collection happened
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptKnown = ADDR_TMP_PERM<T>();

 }

 }

 // return the order
 return ord;
}

static Obj FuncORDER_PERM(Obj self, Obj perm)
{
 RequirePermutation(SELF_NAME, perm);

 if ( TNUM_OBJ(perm) == T_PERM2 ) {
 return ORDER_PERM<UInt2>(perm);
 }
 else {
 return ORDER_PERM<UInt4>(perm);
 }
}

/****************************************************************************
**
*F FuncSIGN_PERM( <self>, <perm> ) . . . . . . . . . . sign of a permutation
**
** 'FuncSIGN_PERM' implements the internal function 'SignPerm'.
**
** 'SignPerm( <perm> )'
**
** 'SignPerm' returns the sign of the permutation <perm>. The sign is +1 if
** <perm> is the product of an *even* number of transpositions, and -1 if
** <perm> is the product of an *odd* number of transpositions. The sign
** is a homomorphism from the symmetric group onto the multiplicative group
** $\{ +1, -1 \}$, the kernel of which is the alternating group.
*/
template <typename T>
static inline Obj SIGN_PERM(Obj perm)
{
 const T * ptPerm; // pointer to the permutation
 Int sign; // sign (result)
 T * ptKnown; // pointer to temporary bag
 UInt len; // length of one cycle
 UInt p, q; // loop variables

 // make sure that the buffer bag is large enough
 UseTmpPerm(SIZE_OBJ(perm));

 // get the pointer to the bags
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptKnown = ADDR_TMP_PERM<T>();

 // clear the buffer bag
 for ( p = 0; p < DEG_PERM<T>(perm); p++ )
 ptKnown[p] = 0;

 // start with sign 1
 sign = 1;

 // loop over all cycles
 for ( p = 0; p < DEG_PERM<T>(perm); p++ ) {

 // if we haven't looked at this cycle so far
 if ( ptKnown[p] == 0 && ptPerm[p] != p ) {

 // find the length of this cycle
 len = 1;
 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
 len++; ptKnown[q] = 1;
 }

 // if the length is even invert the sign
 if ( len % 2 == 0 )
 sign = -sign;

 }

 }

 // return the sign
 return INTOBJ_INT( sign );
}

static Obj FuncSIGN_PERM(Obj self, Obj perm)
{
 RequirePermutation(SELF_NAME, perm);

 if ( TNUM_OBJ(perm) == T_PERM2 ) {
 return SIGN_PERM<UInt2>(perm);
 }
 else {
 return SIGN_PERM<UInt4>(perm);
 }
}

/****************************************************************************
**
*F FuncSMALLEST_GENERATOR_PERM( <self>, <perm> ) . . . . . . . . . . . . . .
*F . . . . . . . smallest generator of cyclic group generated by permutation
**
** 'FuncSMALLEST_GENERATOR_PERM' implements the internal function
** 'SmallestGeneratorPerm'.
**
** 'SmallestGeneratorPerm( <perm> )'
**
** 'SmallestGeneratorPerm' returns the smallest generator of the cyclic
** group generated by the permutation <perm>. That is the result is a
** permutation that generates the same cyclic group as <perm> and is with
** respect to the lexicographical order defined by '\<' the smallest such
** permutation.
*/
template <typename T>
static inline Obj SMALLEST_GENERATOR_PERM(Obj perm)
{
 Obj small; // handle of the smallest gen
 T * ptSmall; // pointer to the smallest gen
 const T * ptPerm; // pointer to the permutation
 T * ptKnown; // pointer to temporary bag
 Obj ord; // order, may be huge
 Obj pow; // power, may also be huge
 UInt len; // length of one cycle
 UInt gcd, s, t; // gcd( len, ord ), temporaries
 UInt min; // minimal element in a cycle
 UInt p, q; // loop variables
 UInt l, n, x, gcd2; // loop variable

 // make sure that the buffer bag is large enough
 UseTmpPerm(SIZE_OBJ(perm));

 // allocate the result bag
 small = NEW_PERM<T>( DEG_PERM<T>(perm) );

 // get the pointer to the bags
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptKnown = ADDR_TMP_PERM<T>();
 ptSmall = ADDR_PERM<T>(small);

 // clear the buffer bag
 for ( p = 0; p < DEG_PERM<T>(perm); p++ )
 ptKnown[p] = 0;

 // we only know that we must raise <perm> to a power = 0 mod 1
 ord = INTOBJ_INT(1); pow = INTOBJ_INT(0);

 // loop over all cycles
 for ( p = 0; p < DEG_PERM<T>(perm); p++ ) {

 // if we haven't looked at this cycle so far
 if ( ptKnown[p] == 0 ) {

 // find the length of this cycle
 len = 1;
 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
 len++; ptKnown[q] = 1;
 }

 // compute the gcd with the previously order ord
 /* Note that since len is single precision, ord % len is to*/
 gcd = len; s = INT_INTOBJ( ModInt( ord, INTOBJ_INT(len) ) );
 while ( s != 0 ) {
 t = s; s = gcd % s; gcd = t;
 }

 // we must raise the cycle into a power = pow mod gcd
 x = INT_INTOBJ( ModInt( pow, INTOBJ_INT( gcd ) ) );

 /* find the smallest element in the cycle at such a positio*/
 min = DEG_PERM<T>(perm)-1;
 n = 0;
 for ( q = p, l = 0; l < len; l++ ) {
 gcd2 = len; s = l;
 while ( s != 0 ) { t = s; s = gcd2 % s; gcd2 = t; }
 if ( l % gcd == x && gcd2 == 1 && q <= min ) {
 min = q;
 n = l;
 }
 q = ptPerm[q];
 }

 // raise the cycle to that power and put it in the result
 ptSmall[p] = min;
 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
 min = ptPerm[min]; ptSmall[q] = min;
 }

 // compute the new order and the new power
 while ( INT_INTOBJ( ModInt( pow, INTOBJ_INT(len) ) ) != n )
 pow = SumInt( pow, ord );
 ord = ProdInt( ord, INTOBJ_INT( len / gcd ) );

 }

 }

 // return the smallest generator
 return small;
}

static Obj FuncSMALLEST_GENERATOR_PERM(Obj self, Obj perm)
{
 RequirePermutation(SELF_NAME, perm);

 if ( TNUM_OBJ(perm) == T_PERM2 ) {
 return SMALLEST_GENERATOR_PERM<UInt2>(perm);
 }
 else {
 return SMALLEST_GENERATOR_PERM<UInt4>(perm);
 }
}

/****************************************************************************
**
*F FuncRESTRICTED_PERM( <self>, <perm>, <dom>, <test> ) . . RestrictedPerm
**
** 'FuncRESTRICTED_PERM' implements the internal function
** 'RESTRICTED_PERM'.
**
** 'RESTRICTED_PERM( <perm>, <dom>, <test> )'
**
** 'RESTRICTED_PERM' returns the restriction of <perm> to <dom>. If <test>
** is set to `true' it is verified that <dom> is the union of cycles of
** <perm>.
*/
template <typename T>
static inline Obj RESTRICTED_PERM(Obj perm, Obj dom, Obj test)
{
 Obj rest;
 T * ptRest;
 const T * ptPerm;
 const Obj * ptDom;
 Int i,inc,len,p,deg;

 // make sure that the buffer bag is large enough
 UseTmpPerm(SIZE_OBJ(perm));

 // allocate the result bag
 deg = DEG_PERM<T>(perm);
 rest = NEW_PERM<T>(deg);

 // get the pointer to the bags
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptRest = ADDR_PERM<T>(rest);

 // create identity everywhere
 for ( p = 0; p < deg; p++ ) {
 ptRest[p]=(T)p;
 }

 if ( ! IS_RANGE(dom) ) {
 if ( ! IS_PLIST( dom ) ) {
 return Fail;
 }
 // domain is list
 ptPerm = CONST_ADDR_PERM<T>(perm);
 ptRest = ADDR_PERM<T>(rest);
 ptDom = CONST_ADDR_OBJ(dom);
 len = LEN_LIST(dom);
 for (i = 1; i <= len; i++) {
 if (IS_POS_INTOBJ(ptDom[i])) {
 p = INT_INTOBJ(ptDom[i]);
 if (p <= deg) {
 p -= 1;
 ptRest[p] = ptPerm[p];
 }
 }
 else {
 return Fail;
 }
 }
 }
 else {
 len = GET_LEN_RANGE(dom);
 p = GET_LOW_RANGE(dom);
 inc = GET_INC_RANGE(dom);
 if (p < 1 || p + inc * (len - 1) < 1) {
 return Fail;
 }
 for (i = p; i != p + inc * len; i += inc) {
 if (i <= deg) {
 ptRest[i - 1] = ptPerm[i - 1];
 }
 }
 }

 if (test==True) {

 T * ptTmp = ADDR_TMP_PERM<T>();

 // cleanout
 for ( p = 0; p < deg; p++ ) {
 ptTmp[p]=0;
 }

 // check whether the result is a permutation
 for (p=0;p<deg;p++) {
 inc=ptRest[p];
 if (ptTmp[inc]==1) return Fail; // point was known
 else ptTmp[inc]=1; // now point is known
 }

 }

 // return the restriction
 return rest;
}

static Obj FuncRESTRICTED_PERM(Obj self, Obj perm, Obj dom, Obj test)
{
 RequirePermutation(SELF_NAME, perm);

 if ( TNUM_OBJ(perm) == T_PERM2 ) {
 return RESTRICTED_PERM<UInt2>(perm, dom, test);
 }
 else {
 return RESTRICTED_PERM<UInt4>(perm, dom, test);
 }
}

/****************************************************************************
**
*F FuncTRIM_PERM( <self>, <perm>, <n> ) . . . . . . . . . trim a permutation
**
** 'TRIM_PERM' trims a permutation to the first <n> points. This can be
## useful to save memory
*/
static Obj FuncTRIM_PERM(Obj self, Obj perm, Obj n)
{
 RequirePermutation(SELF_NAME, perm);
 RequireNonnegativeSmallInt(SELF_NAME, n);
 UInt newDeg = INT_INTOBJ(n);
 UInt oldDeg =
 (TNUM_OBJ(perm) == T_PERM2) ? DEG_PERM2(perm) : DEG_PERM4(perm);
 if (newDeg > oldDeg)
 newDeg = oldDeg;
 TrimPerm(perm, newDeg);
 return 0;
}

/****************************************************************************
**
*F FuncSPLIT_PARTITION( <Ppoints>, <Qnum>,<j>,<g>,<l>)
** <l> is a list [<a>,,<max>] -- needed because of large parameter number
**
** This function is used in the partition backtrack to split a partition.
** The points in the list Ppoints between a (start) and b (end) such
** that Qnum[i^g]=j will be moved at the end.
** At most <max> points will be moved.
** The function returns the start point of the new partition (or -1 if too
** many are moved).
** Ppoints and Qnum must be plain lists of small integers.
*/
template <typename T>
static inline Obj
SPLIT_PARTITION(Obj Ppoints, Obj Qnum, Obj jval, Obj g, Obj lst)
{
 Int a;
 Int b;
 Int max;
 Int blim;
 UInt deg;
 const T * ptPerm;
 Obj tmp;

 a=INT_INTOBJ(ELM_PLIST(lst,1))-1;
 b=INT_INTOBJ(ELM_PLIST(lst,2))+1;
 max=INT_INTOBJ(ELM_PLIST(lst,3));
 blim=b-max-1;

 deg=DEG_PERM<T>(g);
 ptPerm=CONST_ADDR_PERM<T>(g);
 while ( (a<b)) {
 do {
 b--;
 if (b<blim) {
 // too many points got moved out
 return INTOBJ_INT(-1);
 }
 } while (ELM_PLIST(Qnum,
 IMAGE(INT_INTOBJ(ELM_PLIST(Ppoints,b))-1,ptPerm,deg)+1)==jval);
 do {
 a++;
 } while ((a<b)
 &&(!(ELM_PLIST(Qnum,
 IMAGE(INT_INTOBJ(ELM_PLIST(Ppoints,a))-1,ptPerm,deg)+1)==jval)));
 // swap
 if (a<b) {
 tmp=ELM_PLIST(Ppoints,a);
 SET_ELM_PLIST(Ppoints,a,ELM_PLIST(Ppoints,b));
 SET_ELM_PLIST(Ppoints,b,tmp);
 }
 }

 // list is not necc. sorted wrt. \< (any longer)
 RESET_FILT_LIST(Ppoints, FN_IS_SSORT);
 RESET_FILT_LIST(Ppoints, FN_IS_NSORT);

 return INTOBJ_INT(b+1);
}

static Obj
FuncSPLIT_PARTITION(Obj self, Obj Ppoints, Obj Qnum, Obj jval, Obj g, Obj lst)
{
 if (TNUM_OBJ(g)==T_PERM2) {
 return SPLIT_PARTITION<UInt2>(Ppoints, Qnum, jval, g, lst);
 }
 else {
 return SPLIT_PARTITION<UInt4>(Ppoints, Qnum, jval, g, lst);
 }
}

/*****************************************************************************
**
*F FuncDISTANCE_PERMS( <perm1>, <perm2> )
**
** 'DistancePerms' returns the number of points moved by <perm1>/<perm2>
**
*/
template <typename TL, typename TR>
static inline Obj DISTANCE_PERMS(Obj opL, Obj opR)
{
 UInt dist = 0;
 const TL * ptL = CONST_ADDR_PERM<TL>(opL);
 const TR * ptR = CONST_ADDR_PERM<TR>(opR);
 UInt degL = DEG_PERM<TL>(opL);
 UInt degR = DEG_PERM<TR>(opR);
 UInt i;
 if (degL < degR) {
 for (i = 0; i < degL; i++)
 if (ptL[i] != ptR[i])
 dist++;
 for (; i < degR; i++)
 if (ptR[i] != i)
 dist++;
 }
 else {
 for (i = 0; i < degR; i++)
 if (ptL[i] != ptR[i])
 dist++;
 for (; i < degL; i++)
 if (ptL[i] != i)
 dist++;
 }

 return INTOBJ_INT(dist);
}

static Obj FuncDISTANCE_PERMS(Obj self, Obj opL, Obj opR)
{
 UInt type = (TNUM_OBJ(opL) == T_PERM2 ? 20 : 40) + (TNUM_OBJ(opR) == T_PERM2 ? 2 : 4);
 switch (type) {
 case 22:
 return DISTANCE_PERMS<UInt2, UInt2>(opL, opR);
 case 24:
 return DISTANCE_PERMS<UInt2, UInt4>(opL, opR);
 case 42:
 return DISTANCE_PERMS<UInt4, UInt2>(opL, opR);
 case 44:
 return DISTANCE_PERMS<UInt4, UInt4>(opL, opR);
 }
 return Fail;
}

/****************************************************************************
**
*F FuncSMALLEST_IMG_TUP_PERM( <tup>, <perm> )
**
** `SmallestImgTuplePerm' returns the smallest image of the tuple <tup>
** under the permutation <perm>.
*/
template <typename T>
static inline Obj SMALLEST_IMG_TUP_PERM(Obj tup, Obj perm)
{
 UInt res; // handle of the image, result
 const Obj * ptTup; // pointer to the tuple
 const T * ptPrm; // pointer to the permutation
 UInt tmp; // temporary handle
 UInt deg; // degree of the permutation
 UInt i, k; // loop variables

 res = MAX_DEG_PERM4; // ``infty''.

 // get the pointer
 ptTup = CONST_ADDR_OBJ(tup) + LEN_LIST(tup);
 ptPrm = CONST_ADDR_PERM<T>(perm);
 deg = DEG_PERM<T>(perm);

 // loop over the entries of the tuple
 for ( i = LEN_LIST(tup); 1 <= i; i--, ptTup-- ) {
 k = INT_INTOBJ( *ptTup );
 if ( k <= deg )
 tmp = ptPrm[k-1] + 1;
 else
 tmp = k;
 if (tmp<res) res = tmp;
 }

 return INTOBJ_INT(res);

}

static Obj FuncSMALLEST_IMG_TUP_PERM(Obj self, Obj tup, Obj perm)
{
 if ( TNUM_OBJ(perm) == T_PERM2 ) {
 return SMALLEST_IMG_TUP_PERM<UInt2>(tup, perm);
 }
 else {
 return SMALLEST_IMG_TUP_PERM<UInt4>(tup, perm);
 }
}

/****************************************************************************
**
*F OnTuplesPerm( <tup>, <perm> ) . . . . operations on tuples of points
**
** 'OnTuplesPerm' returns the image of the tuple <tup> under the
** permutation <perm>. It is called from 'FuncOnTuples'.
**
** The input <tup> must be a non-empty and dense plain list. This is not
** verified.
*/
template <typename T>
static inline Obj OnTuplesPerm_(Obj tup, Obj perm)
{
 Obj res; // handle of the image, result
 Obj * ptRes; // pointer to the result
 const T * ptPrm; // pointer to the permutation
 Obj tmp; // temporary handle
 UInt deg; // degree of the permutation
 UInt i, k; // loop variables

 // copy the list into a mutable plist, which we will then modify in place
 res = PLAIN_LIST_COPY(tup);
 RESET_FILT_LIST(res, FN_IS_SSORT);
 RESET_FILT_LIST(res, FN_IS_NSORT);
 const UInt len = LEN_PLIST(res);

 // get the pointer
 ptRes = ADDR_OBJ(res) + 1;
 ptPrm = CONST_ADDR_PERM<T>(perm);
 deg = DEG_PERM<T>(perm);

 // loop over the entries of the tuple
 for (i = 1; i <= len; i++, ptRes++) {
 tmp = *ptRes;
 if (IS_POS_INTOBJ(tmp)) {
 k = INT_INTOBJ(tmp);
 if (k <= deg) {
 *ptRes = INTOBJ_INT(ptPrm[k - 1] + 1);
 }
 }
 else if (tmp == NULL) {
 ErrorQuit("OnTuples: must not contain holes", 0, 0);
 }
 else {
 tmp = POW(tmp, perm);
 // POW may trigger a garbage collection, so update pointers
 ptRes = ADDR_OBJ(res) + i;
 ptPrm = CONST_ADDR_PERM<T>(perm);
 *ptRes = tmp;
 CHANGED_BAG( res );
 }
 }

 return res;
}

Obj OnTuplesPerm (
 Obj tup,
 Obj perm )
{
 if ( TNUM_OBJ(perm) == T_PERM2 ) {
 return OnTuplesPerm_<UInt2>(tup, perm);
 }
 else {
 return OnTuplesPerm_<UInt4>(tup, perm);
 }
}

/****************************************************************************
**
*F OnSetsPerm( <set>, <perm> ) . . . . . . . . operations on sets of points
**
** 'OnSetsPerm' returns the image of the tuple <set> under the permutation
** <perm>. It is called from 'FuncOnSets'.
**
** The input <set> must be a non-empty set, i.e., plain, dense and strictly
** sorted. This is not verified.
*/
template <typename T>
static inline Obj OnSetsPerm_(Obj set, Obj perm)
{
 Obj res; // handle of the image, result
 Obj * ptRes; // pointer to the result
 const T * ptPrm; // pointer to the permutation
 Obj tmp; // temporary handle
 UInt deg; // degree of the permutation
 UInt i, k; // loop variables

 // copy the list into a mutable plist, which we will then modify in place
 res = PLAIN_LIST_COPY(set);
 const UInt len = LEN_PLIST(res);

 // get the pointer
 ptRes = ADDR_OBJ(res) + 1;
 ptPrm = CONST_ADDR_PERM<T>(perm);
 deg = DEG_PERM<T>(perm);

 // loop over the entries of the tuple
 BOOL isSmallIntList = TRUE;
 for (i = 1; i <= len; i++, ptRes++) {
--> --------------------

--> maximum size reached

--> --------------------

quality92%

¤ Dauer der Verarbeitung: 0.15 Sekunden (vorverarbeitet) ¤

Versionsinformation zu Columbo

Bemerkung:

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Anfrage:

Dauer der Verarbeitung:

Sekunden

sprechenden Kalenders