/*********************************************************************
*********************************************************************
*********************************************************************
*
*  xperm.c
*
* 	(C) Jose M. Martin-Garcia 2003-2008.
*      jmm@iem.cfmac.csic.es, IEM, CSIC, Madrid, Spain.
*
*	This is free software, distributed under the GNU GPL license.
*      See http://metric.iem.csic.es/Martin-Garcia/xAct/
*
*  These are a collection of C-functions that find Strong Generating
*  Sets, Coset representatives and Double-coset representatives.
*  The algorithms are based on Butler and Portugal et al.
*
*  xperm can be used standalone or from the Mathematica package xPerm.
*
*  20 June - 5 July 2003. Main development.
*   3 September 2004. Eliminated MAX variables.
*   9 November 2005. Corrected treatment of SGS of group D.
*   6 May 2006. All arrays declared dynamically to avoid stack limits.
*     Thanks to Kasper Peeters for spotting and correcting this.
*  25-28 June 2007. Large extension to included multiple dummysets and
*     repeatedsets.
*   2015. Eliminate nested functions. Kasper Peeters.
*
*  Main ideas:
*      - Permutations are represented using Images notation.
*      - Generating sets with m n-permutations are stored as lists of
*	  length m*n.
*	- #define VERBOSE_* to turn on log output.               *PPC*
*
*  Comments:
*	- Permutations are assumed to have degree n>0.
*	- Lists can have length 0 or positive.
*
*  This is ISO C99, not ANSI-C. There are some gcc extensions:
*	- ISO C forbids nested functions
*	- ISO C89 forbids mixed declarations and code
*	- ISO C90 does not support `long long'
*	- ISO C90 forbids variable-size arrays
*
*********************************************************************
*********************************************************************
*********************************************************************/

/* KP */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <vector>
#include "xperm_new.h"
#include <iostream>

/*********************************************************************
*                             PROTOTYPES                            *
*********************************************************************/

/* Output */
void print_perm(int *p, int n, int nl);
void print_array_perm(int *perms, int m, int n, int nl);
void print_list(int *list, int n, int nl);
void print_array(int *array, int m, int n, int nl);

/* Lists */
int equal_list(int *list1, int *list2, int n);
void copy_list(int *list1, int *list2, int n);
int position(int i, int *list, int n);
int position_list(int *matrix, int m, int *row, int n);
void zeros(int *list, int n);
void range(int *list, int n);
void complement(int *all, int al, int *part, int pl, int n,
                int *com, int *cl);
void sort(int *list, int *slist, int l);
void sortB(int *list, int *slist, int l, int *B, int Bl);
int minim(int *list, int n);
int maxim(int *list, int n);
void intersection(int *list1, int l1, int *list2, int l2, int *list,
                  int *l);

/* Permutations */
int isid(int *list, int n );
void product(int *p1, int *p2, int *p, int n);
void inverse(int *p, int *ip, int n);
int onpoints(int point, int *p, int n);
void stable_points(int *p, int n, int *list, int *m);
int first_nonstable_point(int *p, int n);
void nonstable_points(int *list1, int l1, int *GS, int m, int n,
                      int *list2, int *l2);
void stabilizer(int *points, int k, int *GS, int m, int n, int *subGS,
                int *mm);
void one_orbit(int point, int *GS, int m, int n, int *orbit, int *ol);
void all_orbits(int *GS, int m, int n, int *orbits);
void schreier_vector(int point, int *GS, int m, int n, int *nu, int *w);
void trace_schreier(int point, int *nu, int *w, int *perm, int n);

/* Algorithms */
long long int order_of_group(int *base, int bl, int *GS, int m, int n);
int perm_member(int *p, int *base, int bl, int *GS, int m, int n);
void schreier_sims_step(int *base, int bl, int *GS, int m, int n, int i,
                        int *T, int mm, int *newbase, int *nbl, int **newGS, int *nm,
                        int *num);
void schreier_sims(int *base, int bl, int *GS, int m, int n,
                   int *newbase, int *nbl, int **newGS, int *nm, int *num);
void coset_rep(int *p, int n, int *base, int bl, int *GS, int *m,
               int *freeps, int fl, int *cr);
void SGSD(int *vds, int vdsl, int *dummies, int dl, int *mQ,
          int *vrs, int vrsl, int *repes, int rl, int n,
          int firstd, int *KD, int *KDl, int *bD, int *bDl);
void double_coset_rep(int *g, int n, int *base, int bl, int *GS, int m,
                      int *vds, int vdsl, int *dummies, int dl, int *mQ,
                      int *vrs, int vrsl, int *repes, int rl, int *dcr);
void canonical_perm(int *perm,
                    int SGSQ, int *base, int bl, int *GS, int m, int n,
                    int *freeps, int fl, int *dummyps, int dl, int ob, int metricQ,
                    int *cperm);
void canonical_perm_ext(int *perm, int n,
                        int SGSQ, int *base, int bl, int *GS, int m,
                        int *frees, int fl,
                        int *vds, int vdsl, int *dummies, int dl, int *mQ,
                        int *vrs, int vrsl, int *repes, int rl,
                        int *cperm);


/*********************************************************************
*                         PRINTING FUNCTIONS                        *
*********************************************************************/

/**********************************************************************/

/* print_perm. JMM, 22 June 2003
*
* This function prints a permutation p of degree n. If nl=1(0) adds
* (does not) a newline. */

void print_perm(int *p, int n, int nl)
	{
	int i;
	if (isid(p,n)) printf("id");
	else {
		printf("(");
		printf("%d", p[0]);          /* No comma */
		for (i=1; i<n; i++) {
			printf(",%d", p[i]); /* Comma separated */
			}
		printf(")");
		}
	if (nl) printf("\n");
	}

/**********************************************************************/

/* print_array_perm. JMM, 22 June 2003
*
* This function prints an array of m n-permutations.
* If nl=1 (0) adds (does not) a newline after each row.
* There are no commas between permutations. */

void print_array_perm(int *perms, int m, int n, int nl)
	{
	int j;
	printf("{");
	if (nl) printf("\n");
	for(j=0; j<m; j++) {
		printf(" ");
		print_perm(perms+j*n, n, nl);
		}
	if (nl) printf("}\n");
	else printf(" }\n");
	}

/**********************************************************************/

/* print_list. JMM, 22 June 2003
*
* This function prints a list of length n in curly brackets.
* If nl=1 (0) it adds (does not) a newline. */

void print_list(int *list, int n, int nl)
	{
	int i;
	printf("{");
	if (n>0) printf("%d", list[0]);            /* No comma */
	for (i=1; i<n; i++) printf(",%d", list[i]);/* Comma separated */
	printf("}");
	if (nl) printf("\n");
	}

/**********************************************************************/

/* print_array. JMM, 22 June 2003
*
* This function prints an array of dimensions m x n in curly brackets.
* If nl=1 (0) adds (does not) a newline after each row.
* There are no commas between lists. */

void print_array(int *array, int m, int n, int nl)
	{
	int j;
	printf("{");
	if(nl) printf("\n");
	for(j=0; j<m; j++) {
		printf(" ");
		print_list(array+j*n, n, nl);
		}
	if (!nl) printf(" ");
	printf("}\n");
	}

/**********************************************************************/


/*********************************************************************
*                          GENERIC FUNCTIONS                        *
*********************************************************************/

/**********************************************************************/

/* equal_list. JMM, 27 June 2003
*
* Checks that list1 and list2 (both with length n) are equal.
* If n=0 the result is 1.
*/

/* Old code
int equal_list(int *list1, int *list2, int n) {
	while(n--) { * Run from n-1 to 0 *
	if(*(list1+n) != *(list2+n)) return(0); * different *
			}
	return(1); * equal *
			}
*/

/* KP, 7 May 2006 */
int equal_list(int *list1, int *list2, int n)
	{
	if (n==0) return 1;
	if (memcmp(list1, list2, n*sizeof(int))==0) return 1;
	else return 0;
	}

/**********************************************************************/

/* copy_list. JMM, 27 June 2003
*
* Copies first n elements of list1 onto list2. If n=0 nothing is done.
*/

/* Old code
void copy_list(int *list1, int *list2, int n) {
	while(n--) *(list2+n) = *(list1+n);
			}
*/

/* KP, 7 May 2006 */
/* KP, 26 Nov 2007: memcpy changed to memmove */
void copy_list(int *list1, int *list2, int n)
	{
	if (n==0) return;
	memmove(list2, list1, n*sizeof(int));
	}

/**********************************************************************/

/* position. JMM, 20 June 2003
*
* This function finds the position of the integer i in list (length n),
* counting from 1 to n, or else gives 0 if i is not in list. Note that
* it starts from the end so that if i appears several times in list it
* will only find the last one. If n=0 the result is 0.
*
* More than 60% of time in the Schreier-Sims algorithm is spent here!
*
* Note: it is not possible to use  memchr  because it only works with
* characters (integers up to 255).
*/

int position(int i, int *list, int n)
	{
	while(n--) {
		if (*(list+n) == i) return(n+1);
		}
	return(0);
	}

/**********************************************************************/

/* position_list. JMM, 28 June 2003
*
* This function gives the position of list row (length n) in the
* matrix of dimensions m, n. If row is not present in matrix, 0
* is returned.
*/

int position_list(int *matrix, int m, int *row, int n)
	{
	while(m--) {
		if(equal_list(matrix+m*n, row, n)) return(m+1);
		}
	return(0);
	}

/**********************************************************************/

/* zeros, range. JMM, 22 June 2003
*
* These functions generate lists of n zeros or the range 1...n,
* respectively, storing the result in list.
*/

/* Old code
void zeros(int *list, int n) {
	while(n--) *(list+n) = 0;
			}
*/

/* KP, 7 May 2006 */
void zeros(int *list, int n)
	{
	if (n==0) return;
	memset(list, 0, n*sizeof(int));
	}

void range(int *list, int n)
	{
	while(n--) *(list+n) = n+1;
	}

/**********************************************************************/

/* complement. JMM, 28 June 2003
*
* This function gives all sublists of length n in list all (length al)
* that are not members of list part (length pl), storing the result
* in list com (length cl).
* We assume that enough space (typically al n-lists) has been already
* allocated for com.
*/

void complement(int *all, int al, int *part, int pl, int n,
                int *com, int *cl)
	{
	int i;
	*cl = 0;
	for(i=0; i<al; i++) {
		if (!position_list(part, pl, all+i*n, n)) {
			copy_list(all+i*n, com+(*cl)*n, n);
			(*cl)++;
			}
		}
	}

/**********************************************************************/

/* sort. JMM, 30 June 2003
*
* Sorts list (length l) storing the result in slist (length l, too).
* There are l*(l-1)/2 element comparisons. This could be improved.
*/

void sort(int *list, int *slist, int l)
	{

	int i,j;
	int tmp, mini; /* mini is a position, not an element of list */

#ifdef VERBOSE_LISTS						/*PPC*/
	printf("sort: Sorting list ");
	print_list(list, l, 1);	/*PPC*/
#endif								/*PPC*/
	copy_list(list, slist, l);
	for (i=0; i<l-1; i++) {
		/* Assume points at positions 0...i-1 are already sorted */
		/* Compare point at position i with points to the right */
		mini = i;
		for(j=i+1; j<l; j++) {
			if ( slist[j] < slist[mini] ) mini = j;
			}
		/* Swap points i and mini */
		tmp = slist[i];
		slist[i] = slist[mini];
		slist[mini] = tmp;
		}
#ifdef VERBOSE_LISTS						/*PPC*/
	printf("sort: with result ");
	print_list(slist, l, 1);	/*PPC*/
#endif								/*PPC*/
	}

/* sortB. JMM, 30 June 2003
*
* Sorts list (length l) according to the order given by B (length Bl),
* storing the result in slist (length l, too). Elements not in B
* are pushed to the end, and sorted using sort.
*/

void sortB(int *list, int *slist, int l, int *B, int Bl)
	{

	int sl;
	int *tmp=  (int*)malloc(l*sizeof(int)), tmpl;
	int *stmp= (int*)malloc(l*sizeof(int));

#ifdef VERBOSE_LISTS						/*PPC*/
	printf("sortB: Sorting list ");
	print_list(list, l, 1);	/*PPC*/
	printf("sortB: using base ");
	print_list(B, Bl, 1);	/*PPC*/
#endif								/*PPC*/
	/* Elements of list in B, keeping order of B, moved to slist */
	intersection(B, Bl, list, l, slist, &sl);
	/* Other elements of list are appended to slist */
	complement(list, l, B, Bl, 1, tmp, &tmpl);
	/* Check that tmpl+sl==l */
	if(tmpl+sl != l) printf("Error in sortB\n");
	/* Sort the latter integers using sort */
	sort(tmp, stmp, tmpl);
	copy_list(stmp, slist+sl, tmpl);
#ifdef VERBOSE_LISTS						/*PPC*/
	printf("sortB: with result ");
	print_list(slist, l, 1);	/*PPC*/
#endif								/*PPC*/
	/* Free allocated memory */
	free(tmp);
	free(stmp);

	}

/**********************************************************************/

/* minim, maxim. JMM, 2 July 2003
*
* These functions find the minimum and maximum elements of a given
* list of length n, using normal integer ordering. Names min, max are
* reserved.
*/

int minim(int *list, int n)
	{

	int m=list[n-1];

	while(n--) {
		if (list[n]<m) m=list[n];
		}
	return(m);

	}

int maxim(int *list, int n)
	{

	int m=list[n-1];

	while(n--) {
		if (list[n]>m) m=list[n];
		}
	return(m);

	}

/**********************************************************************/

/* intersection. JMM, 3 July 2003
*
* This function returns in list (length l) the intersection of the
* elements of lists list1 (length l1) and list2 (length l2). Note that
* we return unique (that is, non-repeated) elements. Returned elements
* keep the original order in list1.
* We assume that enough space (typical l1 or l2 integers) has been
* already allocated for list.
*/

void intersection(int *list1, int l1, int *list2, int l2, int *list,
                  int *l)
	{

	int i, j, a;

	*l = 0;
	for(i=0; i<l1; i++) { /* Range over elements of list1 */
		a = list1[i];
		for(j=0; j<l2; j++) { /* Range over elements of list2 */
			/* a is in list2 and not yet in list, append */
			if ((list2[j] == a) && (!position(a, list, *l)))
				list[(*l)++] = a;
			}
		}

	}

/**********************************************************************/


/*********************************************************************
*                        PERMUTATIONS FUNCTIONS                     *
*********************************************************************/

/**********************************************************************/

/* isid. JMM, 20 June 2003
*
* This function detects whether the n-permutation p is the identity
* (returning 1) or not (returning 0). */

int isid(int *p, int n )
	{

	while(n--) {
		if (*(p+n)!=n+1) return(0); /* Not the identity */
		}
	return(1); /* Identity */

	}

/**********************************************************************/

/* product. JMM, 20 June 2003
*
* Product of n-permutations p1 and p2 (from left to right) storing the
* result in the n-permutation p. Note that the points of the
* permutations go from 1 to n, while we want to add from 0 to n-1 to
* the pointer p2. That's why we need to add -1.
* This function is highly efficient, but rather criptic. */

void product(int *p1, int *p2, int *p, int n)
	{

	while(n--) *(p++) = *(p2-1+*(p1++));

	/* Example:
	* Suppose we have p1=(3,2,1) and p2=(3,1,2) of degree n=3.
	* The following operations take place, in this order:
	*   n=3, n--
	*   *(p1++)=3
	*   *(p2-1+3)=*(p2+2)=2
	*   *(p++)=2  Computed first element of p
	*   n=2, n--
	*   *(p1++)=2
	*   *(p2-1+2)=*(p2+1)=1
	*   *(p++)=1  Computed second element of p
	*   n=1, n--
	*   *(p1++)=1
	*   *(p2-1+1)=*(p2)=3
	*   *(p++)=3  Computed third element of p
	*   n=0, END
	* Therefore the final permutation returned is p=(2,1,3).
	*/

	}

/* Taken out KP */
/* TAB1 is an element of S; TAB2 is an element of D */
void F2(int *TAB1, int *g, int *TAB2, int *sgd, int n)
	{
	int *tmp= (int*)malloc(n*sizeof(int));

	product(TAB1, g, tmp, n);
	product(tmp, TAB2, sgd, n);

	free(tmp);
	} /* End of function F2 */

/**********************************************************************/

/* inverse. JMM, 20 June 2003
*
* Inverse of the n-permutation p, storing the result in the
* n-permutation ip. Again we need to add -1 to the pointer ip.
* Conversely, note that we need to add +1 to n because n ranges from
* 0 to n-1 in the while loop. */

void inverse(int *p, int *ip, int n)
	{

	while(n--) *(ip-1+*(p+n)) = n+1 ;

	}

/**********************************************************************/

/* onpoints. JMM, 20 June 2003
*
* Image of point under the n-permutation p. If point is larger than
* the degree n we return point. */

int onpoints(int point, int *p, int n)
	{

	if (point <= n) return(*(p-1+point));
	else return(point);

	}

/**********************************************************************/

/* stable_points. JMM, 20 June 2003
*
* This function finds the points which remain stable under the
* n-permutation p, storing the result in list of length m (<=n).
* With n=0 we get an empty list (m=0).
* We assume that enough space (typically n integers) has been already
* allocated for list.
*/

void stable_points(int *p, int n, int *list, int *m)
	{

	int i;

	*m = 0;
	for(i=1; i<=n; i++) {
		if ( onpoints(i, p, n) == i ) list[(*m)++] = i;
		}

	}

/**********************************************************************/

/* first_nonstable_point. JMM, 30 June 2003
*
* This function gives the smallest moved point by the permutation p.
* If no point is moved the result is 0. If n=0 the result is 0.
* This is a restricted form of stable_points.
*/

int first_nonstable_point(int *p, int n)
	{

	int i;

	for(i=1; i<=n; i++) {
		if (onpoints(i, p, n) != i) return(i);
		}
	return(0);
	}

/**********************************************************************/

/* nonstable_points. JMM, 30 June 2003
*
* This function gives a set of points (list2, length l2) such that
* none of the permutations in GS fixes them all. The first l1 points
* of list2 will be those of list1, even if they are all stable under
* GS. If m=0 the resulting list2 is just list1.
* We assume that enough space (typically n integers) has been already
* allocated for list2.
*/

void nonstable_points(int *list1, int l1, int *GS, int m, int n,
                      int *list2, int *l2)
	{

	int i, j;

	copy_list(list1, list2, l1);
	*l2 = l1; /* Initialize list2 */
	for(j=0; j<m; j++) { /* Range over permutations of GS */
		int stable=1;
		for(i=0; i<*l2; i++) {
			if (onpoints(list2[i], GS+j*n, n) != list2[i]) {
				stable=0;
				break;
				}
			}
		/* If all points already in list2 are stable under the
			permutation, append the smallest nonstable point */
		if (stable) {
			list2[*l2] = first_nonstable_point(GS+j*n, n);
			(*l2)++;
			}
		}
	}

/**********************************************************************/

/* stabilizer. JMM, 22 June 2003
*
* This function finds the permutations that stabilize the set of k
* points among the GeneratingSet GS with m n-permutations. With m=0
* we return an empty subset (mm=0). Recall that the stabilizer of a
* generating set of a group is not, in general, a generating set for
* the corresponding stabilizer of that group.
* We assume that enough space (typically m*n integers) has been already
* allocated for subGS.
*/

void stabilizer(int *points, int k, int *GS, int m, int n,
                int *subGS, int *mm)
	{

	int i, j;
	int stable;

	*mm=0;
	for(j=0; j<m; j++) {
		stable=1;
		for(i=0; i<k; i++) {
			if (onpoints(points[i],GS+n*j,n) != points[i]) {
				stable = 0;
				break;
				}
			}
		if (stable) {
			copy_list(GS+n*j, subGS+(*mm)*n, n);
			++(*mm);
			}
		}

	}

/**********************************************************************/

/* one_orbit. JMM, 22 June 2003
*
* Orbit of a given point under a GeneratingSet GS with m
* n-permutations. The result is stored in orbit (length ol).
* With m=0 the orbit is just {point}.
* We assume that enough space (typically n integers) has been already
* allocated for orbit.
*/

void one_orbit(int point, int *GS, int m, int n, int *orbit, int *ol)
	{

	int np=0;  /* Index of current element in the orbit */
	int gamma; /* Current element in the orbit */
	int mp;    /* Index of current permutation in GS */
	int newgamma;

	orbit[0] = point;
	*ol = 1;
	while(np < *ol) {
		gamma = orbit[np];
		for(mp=0; mp<m; mp++) {
			newgamma = onpoints(gamma, GS+mp*n, n);
			if (!position(newgamma, orbit, *ol))
				orbit[(*ol)++] = newgamma;
			}
		np++;
		}

	}

/**********************************************************************/

/* all_orbits. JMM, 1 July 2003
*
* This function returns all orbits of the generating set GS (with
* m n-permutations). Note that we use a very special notation: the
* result is a vector orbits of length n such that all points marked
* with 1 belong to the same orbit; all points marked with 2 belong to
* another orbit, and so on. With m=0 we return the list {1,2,...,n}.
*/

void all_orbits(int *GS, int m, int n, int *orbits)
	{

	int i, j;                                   /* Counters */
	int *orbit= (int*)malloc(n*sizeof(int)), ol;/* Computed orbit */
	int orbit_index=1;                                   /* Orbit index */

	/* Initialize orbits */
	memset(orbits, 0, n*sizeof(int));

	/* Compute orbits */
	for (i=1; i<=n; i++) { /* Points */
		if(orbits[i-1]==0) { /* New orbit condition */
			/* Compute new orbit */
			one_orbit(i, GS, m, n, orbit, &ol);
			/* Mark points of new orbit with index or */
			for (j=0; j<ol; j++) orbits[orbit[j]-1] = orbit_index;
			/* Increment index for next orbit */
			orbit_index++;
			}
		}

	/* Free allocated memory */
	free(orbit);

	}

/**********************************************************************/

/* one_schreier_orbit. JMM, 22 June 2003
*
* Orbit of point under the GeneratingSet GS with m n-permutations.
* The result is stored in orbit and the vectors nu of permutations and
* w of backward points. If init=1 both nu and w are reset to 0.
*
* Note: the profiler shows that roughly half of the time in this
* function is spent in the subroutine `position'.
*/

void one_schreier_orbit(int point, int *GS, int m, int n,
                        int *orbit, int *ol, int *nu, int *w, int init)
	{

	int np;    /* Index of current element in the orbit */
	int gamma; /* Current element in the orbit */
	int mp;    /* Index of current permutation in GS */
	int *perm= (int*)malloc(n*sizeof(int));
	int newgamma;

	/* Initialize schreier with zeros if required */
	memset(orbit, 0, n*sizeof(int));
	if (init) {
		memset(nu, 0, n*n*sizeof(int));
		memset(w, 0, n*sizeof(int));
		}
	/* First element of orbit. There is no backward pointer */
	orbit[0] = point;
	*ol = 1;
	/* Other elements of orbit */
	np = 0;
	while(np < *ol) {
		gamma = orbit[np];
		for(mp=0; mp<m; mp++) {
			copy_list(GS+mp*n, perm, n);
			newgamma = onpoints(gamma, perm, n);
			if (position(newgamma, orbit, *ol));
			else {
				/* Append to orbit */
				orbit[(*ol)++] = newgamma;
				/* Perm moving gamma to newgamma */
				copy_list(perm, nu+(newgamma-1)*n, n);
				/* Gamma backward pointer of newgamma */
				*(w+newgamma-1) = gamma;
				}
			}
		np++;
		}
	/* Free allocated memory */
	free(perm);

	}

/**********************************************************************/

/* schreier_vector. JMM, 22 June 2003
*
* All orbits of a given GeneratingSet GS with m n-permutations, with
* combined vectors of backward permutations and pointers. Note that we
* do not return the orbits information, though it can be reconstructed
* from w and nu.
*/

void schreier_vector(int point, int *GS, int m, int n, int *nu, int *w)
	{

	int i;    /* Point counter (from 1 to n) */
	int *orbit=      (int*)malloc(n*sizeof(int));
	int *usedpoints= (int*)malloc(n*sizeof(int));
	int j=0;  /* Counter of used points */
	int ol;

	/* First orbit */
	one_schreier_orbit(point, GS, m, n, orbit, &ol, nu, w, 1);
	while(ol--) usedpoints[j++] = orbit[ol];
	/* Other orbits */
	for(i=1; i<=n; i++) {
		if (!position(i, usedpoints, j)) {
			one_schreier_orbit(i, GS, m, n, orbit, &ol, nu, w, 0);
			while(ol--) usedpoints[j++] = orbit[ol];
			}
		}
	/* Free allocated memory */
	free(orbit);
	free(usedpoints);

	}

/**********************************************************************/

/* trace_schreier. JMM, 22 June 2003
*
* This function traces the Schreier vector (orbits, nu, w) with the
* point given, returning in perm a permutation wich moves the first
* point of the corresponding orbit to point.
*/

void trace_schreier(int point, int *nu, int *w, int *perm, int n)
	{

	int i;
	int found=0;
	int *newperm= (int*)malloc(n*sizeof(int));

	for(i=0; i<n; i++) {
		if (*(w+point-1) == 0) {
			range(perm, n);
			found=1;
			break;
			}
		}
	if (!found) {
		trace_schreier(*(w+point-1), nu, w, newperm, n);
		product(newperm, nu+(point-1)*n, perm, n);
		}
	free(newperm);

	}

/**********************************************************************/

/* order_of_group. JMM, 27 June 2003
*
* This function calculates the order of the group generated by the SGS
* given by base (with bl points) and GS (with m n-permutations).
*
* JMM, 12 March 2006: changed type to long long int (8 bytes).
* This allows numbers up to ~ 10^19.
*/

long long int order_of_group(int *base, int bl, int *GS, int m, int n)
	{

	if (m==0) return(1);
	else {
		int *stab=  (int*)malloc(m*n*sizeof(int)), sl;
		int *orbit= (int*)malloc(  n*sizeof(int)), ol;
		one_orbit(base[0], GS, m, n, orbit, &ol);
		stabilizer(base, 1, GS, m, n, stab, &sl);
		long long int ret=ol* order_of_group(base+1,bl-1,stab,sl,n);
		free(stab);
		free(orbit);
		return ret;
		}
	}

/**********************************************************************/

/* perm_member. JMM, 27 June 2003
*
* Checks whether permutation p belongs (returns 1) or not (returns 0)
* to the group generated by the SGS given by base (length bl) and GS
* (m n-permutations). If bl=0 or m=0 the result is  isid(p, n) .
*
* JMM, 26 June 2007: modified behaviour of m==0. It gave ret=1 before.
*/

int perm_member(int *p, int *base, int bl, int *GS, int m, int n)
	{

	if (bl==0 || m==0) return( isid(p, n) );
	else {
		int *pp=    (int*)malloc(  n*sizeof(int));
		int *ip=    (int*)malloc(  n*sizeof(int));
		int *orbit= (int*)malloc(  n*sizeof(int)), ol;
		int *w=     (int*)malloc(  n*sizeof(int));
		int *nu=    (int*)malloc(n*n*sizeof(int));
		int *stab=  (int*)malloc(m*n*sizeof(int)), sl;
		int point, ret;

		one_schreier_orbit(base[0], GS,m,n, orbit,&ol, nu,w, 1);
		point = onpoints(base[0], p, n);
		if (position(point, orbit, ol)) {
			trace_schreier(point, nu, w, pp, n);
			inverse(pp, ip, n);
			product(p, ip, pp, n);
			stabilizer(base, 1, GS, m, n, stab, &sl);
			ret = perm_member(pp, base+1, bl-1, stab,sl,n);
			}
		else ret = 0;

		free(pp);
		free(ip);
		free(orbit);
		free(w);
		free(nu);
		free(stab);

		return(ret);
		}

	}

/**********************************************************************/

/* schreier_sims. JMM, 28-30 June 2003
*
* This function computes a Strong Generating Set for the group of
* permutations generated by GS (with m n-permutations). It is
* possible to give the first bl (possibly 0) elements of the base,
* and the code computes all the others. There are possibly (and
* probably) redundant points in the final base. n=0 is not consistent
* here.
* We assume that enough space (at least, and typically, m*n integers)
* has been already allocated for newGS. We also assume that newGS can
* be reallocated. That's why we do not send a pointer, but a pointer
* to that pointer. That is, it needs to be reallocated in a different
* subroutine, and that cannot be done with a normal pointer!
*/

void schreier_sims(int *base, int bl, int *GS, int m, int n,
                   int *newbase, int *nbl, int **newGS, int *nm, int *num)
	{

#ifdef VERBOSE_SCHREIER						/*PPC*/
	printf("******** SCHREIER-SIMS ALGORITHM ********\n");	/*PPC*/
#endif								/*PPC*/
	/* Note the cycle: base -> newbase -> base2 -> newbase -> ... */

	/* Copy base into newbase, adding more points if needed */
	nonstable_points(base, bl, GS, m, n, newbase, nbl);
#ifdef VERBOSE_SCHREIER						/*PPC*/
	printf("Original base:");
	print_list(base, bl, 1);	/*PPC*/
	printf("New base:");
	print_list(newbase, *nbl, 1);	/*PPC*/
#endif								/*PPC*/
	/* Initialize newGS=GS */
	copy_list(GS, *newGS, m*n);
	*nm = m;
	if (*nbl==0) { /* Problem. Return input sets  */
		copy_list(base, newbase, bl);
		*nbl = bl;
		return;
		}

	/* Allocate memory for intermediate SGS and stabilizer */
	int i;                                         /* Base index */
	int *base2= (int*)malloc(  n*sizeof(int)), bl2;/* Interm base */
	int *GS2=   (int*)malloc(m*n*sizeof(int)), m2; /* Interm GS */
	int *stab=  (int*)malloc(m*n*sizeof(int)), mm; /* Stabilizer */

	/* Main loop */
	m2 = *nm; /* Initially GS2 will be just newGS */
	for(i=(*nbl); i>0; i--) {
#ifdef VERBOSE_SCHREIER						/*PPC*/
		printf("\nComputing SGS for H^(%d)\n", i-1);		/*PPC*/
#endif								/*PPC*/
		if (*nm > m2) { /* Reallocate GS2 and stab */
			GS2 =  (int*)realloc(GS2,  (*nm)*n*sizeof(int));
			stab = (int*)realloc(stab, (*nm)*n*sizeof(int));
			}
		/* Copy newbase into base2 */
		copy_list(newbase, base2, *nbl);
		bl2=*nbl;
		copy_list(*newGS, GS2, (*nm)*n);
		m2=*nm;
		/* Compute newbase from base2 */
		stabilizer(base2, i-1, GS2, m2, n, stab, &mm);
		schreier_sims_step(base2, bl2, GS2, m2, n, i,
		                   stab, mm, newbase, nbl, newGS, nm, num);
		}

	/* Free allocated memory */
	free(base2);
	free(GS2);
	free(stab);

#ifdef VERBOSE_SCHREIER						/*PPC*/
	printf("************ END OF ALGORITHM ***********\n");	/*PPC*/
#endif								/*PPC*/
	}

/* Intermediate function for schreier_sims.
*
* Assuming that S^(i-1) is a GenSet for H^(i-1) and that [B, S^(i)]
* are a StrongGenSet for H^(i), find a StrongGenSet for H^(i-1).
* T is the set of additional generators in S^(i-1) since the previous
* call to the procedure with the present value of i. Assume that a
* StrongGenSet of < S^(i-1) - T > (the previous value of H^(i-1)) is
* included in B and S. The present value of nu^(i-1) must be an
* extension of the previous value.
*
* base (length bl): tentative base for the whole group
* GS (m n-permutations): tentative StrongGenSet for the whole group
* i: current order in the hierarchy of stabilizers
*
* Sizes on input:
*	i <= bl <= n   Not changed
*      mm <= m        Not changed
*      nbl = bl       nbl will increase
*      nm = m         nm will increase
*/

void schreier_sims_step(int *base, int bl, int *GS, int m, int n,
                        int i, int *T, int mm,
                        int *newbase, int *nbl, int **newGS, int *nm, int *num)
	{

	/* Declarations */
	/* Counters */
	int c, j=0, jj, level;
	/* Intermediate permutations */
	int *p=   (int*)malloc(n*sizeof(int));
	int *ip=  (int*)malloc(n*sizeof(int));
	int *pp=  (int*)malloc(n*sizeof(int));
	int *ppp= (int*)malloc(n*sizeof(int));
	/* Stabilizer of base[1...i-1] */
	int *Si= (int*)malloc(m*n*sizeof(int)), Sil;
	/* Old stabilizer. Here we could use mm*n rather than m*n */
	int *oldSi= (int *)malloc(m*n*sizeof(int)), oldSil;
	/* Orbit of base[i] */
	int *Deltai= (int*)malloc(  n*sizeof(int)), Deltail;
	int *w=      (int*)malloc(  n*sizeof(int));
	int *nu=     (int*)malloc(n*n*sizeof(int));
	/* Old orbit */
	int *oldDeltai= (int*)malloc(  n*sizeof(int)), oldDeltail;
	int *oldw=      (int*)malloc(  n*sizeof(int));
	int *oldnu=     (int*)malloc(n*n*sizeof(int));
	/* Generators to check */
	int *genset= (int*)malloc(m*n*sizeof(int)), gensetl;
	/* Loops */
	int gamma, gn, sn;
	int *s= (int*)malloc(n*sizeof(int));
	int *g= (int*)malloc(n*sizeof(int));
	/* Stabilizer */
	int *stab =  (int*)malloc(m*n*sizeof(int)), stabl;
	int *stabps= (int*)malloc(  n*sizeof(int)), stabpsl;

#ifdef VERBOSE_SCHREIER						/*PPC*/
	printf("******** schreier_sims_step ********\n");	/*PPC*/
	printf("base:");
	print_list(base, bl, 1);		/*PPC*/
	printf("GS (%d perms of degree %d):", m, n);		/*PPC*/
	print_array_perm(GS, m, n, 1);				/*PPC*/
#endif								/*PPC*/
	/* Initialize newbase=base and newGS=GS (and their lengths) */
	/* They are already equal on input. Always true? */
	copy_list(base, newbase, bl);
	*nbl = bl;
	copy_list(GS, *newGS, m*n);
	*nm = m;
	/* Original generating sets. We get Sil<=m and oldSil<=mm */
	stabilizer(base, i-1, GS, m, n, Si, &Sil);
#ifdef VERBOSE_SCHREIER						/*PPC*/
	printf("Stabilizer of first %d points of base ", i-1);	/*PPC*/
	print_list(base, i-1, 0);
	printf(" :\n");		/*PPC*/
	print_array_perm(Si, Sil, n, 1);			/*PPC*/
#endif								/*PPC*/
	complement(Si, Sil, T, mm, n, oldSi, &oldSil);
#ifdef VERBOSE_SCHREIER						/*PPC*/
	printf("Previous stabilizer of %d points:\n", i-1);	/*PPC*/
	print_array_perm(oldSi, oldSil, n, 1);			/*PPC*/
#endif								/*PPC*/
	/* Basic orbits */
	one_schreier_orbit(base[i-1], Si, Sil, n,
	                   Deltai, &Deltail, nu, w, 1);
	one_schreier_orbit(base[i-1], oldSi, oldSil, n,
	                   oldDeltai, &oldDeltail, oldnu, oldw, 1);
	/* Check that Deltai is an extension of oldDeltai */
	for(c=0; c<n; c++) {
		if(w[c]!=oldw[c] && oldw[c]!=0) {
#ifdef VERBOSE_SCHREIER						/*PPC*/
			printf("Deltai[%d] modified to match oldDeltai\n", c);  /*PPC*/
#endif								/*PPC*/
			copy_list(oldnu+c*n, nu+c*n, n);
			w[c] = oldw[c];
			}
		}
#ifdef VERBOSE_SCHREIER						/*PPC*/
	printf("Orbit: ");
	print_list(Deltai, Deltail, 1);	/*PPC*/
#endif								/*PPC*/

	/* Loop gn over elements gamma of basic orbit */
	for(gn=0; gn<Deltail; gn++) {

		gamma = Deltai[gn];

#ifdef VERBOSE_SCHREIER						/*PPC*/
		printf("  gamma=%d\n", gamma);				/*PPC*/
#endif								/*PPC*/
		/* In both cases here we get gensetl<=m */
		if (position(gamma, oldDeltai, oldDeltail)) {
			copy_list(T, genset, mm*n);
			gensetl=mm;
			}
		else {
			copy_list(Si, genset, Sil*n);
			gensetl=Sil;
			}

		/* Loop sn over generators s in genset */
		for (sn=0; sn<gensetl; sn++) {

			copy_list(genset+sn*n, s, n);
			(*num)++;

			/* Compute Schreier generator */
			trace_schreier(gamma, nu,w, p,n);
			trace_schreier(onpoints(gamma,s,n), nu,w, pp,n);
			inverse(pp, ip, n);
			product(p, s, ppp, n);
			product(ppp, ip, g, n);
#ifdef VERBOSE_SCHREIER						/*PPC*/
			printf("    g(%d)=", *num);
			print_perm(g, n, 1);	/*PPC*/
#endif								/*PPC*/

			/* Compute stabilizer. Reallocate to maximum size */
			stab = (int*)realloc(stab, (*nm)*n*sizeof(int));
			stabilizer(newbase, i, *newGS, *nm, n, stab, &stabl);
			/* If g is not in subgroup H^(i) */
			if(!isid(g, n)) {
				if ((stabl==0)||(!perm_member(g, newbase+i,*nbl-i, stab,stabl, n))) {
#ifdef VERBOSE_SCHREIER						/*PPC*/
					printf("      g not in H^(%d)\n", i);			/*PPC*/
#endif								/*PPC*/
					/* Enlarge newGS. We need reallocation */
					*newGS = (int*)realloc(*newGS, ((*nm)+1)*n*sizeof(int));
					copy_list(g, (*newGS)+(*nm)*n, n);
					(*nm)++;
					/* Extend newbase if needed, so that no strong
					* generator fixes all base points. */
					stable_points(g, n, stabps, &stabpsl);
#ifdef VERBOSE_SCHREIER						/*PPC*/
					printf("      Stable under g: "); 			/*PPC*/
					print_list(stabps, stabpsl, 1);				/*PPC*/
#endif								/*PPC*/
					/* If g moves a point of newbase then set j */
					for(jj=0; jj<*nbl; jj++) {
						if(!position(newbase[jj], stabps, stabpsl)) {
							j = jj+1;
#ifdef VERBOSE_SCHREIER						/*PPC*/
							printf("      g moves %d in newbase\n", newbase[jj]);	/*PPC*/
#endif								/*PPC*/
							break;
							}
						/* else choose a new point */
						j = *nbl+1;
						}
					if (j == *nbl+1) {
						for(jj=1; jj<=n; jj++) {
							if(!position(jj, stabps, stabpsl) &&
							      !position(jj, newbase, *nbl)) {
								newbase[*nbl]=jj;
								(*nbl)++;
#ifdef VERBOSE_SCHREIER						/*PPC*/
								printf("      Point %d added to newbase\n", jj);	/*PPC*/
#endif								/*PPC*/
								break;
								}
							}
						}
					/* Ensure we still have a base and SGS for
					* H^(i+1) */
					for(level=j; level>i; level--) {
#ifdef VERBOSE_SCHREIER						/*PPC*/
						printf("\nEnsuring H^(%d) at level %d\n", i+1, level);	/*PPC*/
#endif								/*PPC*/
						schreier_sims_step(newbase,*nbl,*newGS,*nm, n,
						                   level, g, 1,
						                   newbase,nbl,newGS,nm,num);
						}
#ifdef VERBOSE_SCHREIER						/*PPC*/
					printf("***** Finished check of H(%d) ******\n\n", i+1);/*PPC*/
#endif								/*PPC*/
					}
				}
			}
		}

	/* Free allocated memory */
	free(p);
	free(ip);
	free(pp);
	free(ppp);
	free(Si);
	free(oldSi);
	free(Deltai);
	free(w);
	free(nu);
	free(oldDeltai);
	free(oldw);
	free(oldnu);
	free(genset);
	free(s);
	free(g);
	free(stab);
	free(stabps);

	}

/**********************************************************************
*
* Notations and comments:
*  - In xTensor, a permutation g corresponding to a given tensor
*    configuration is understood as acting on index-numbers giving
*    tensor-slot numbers. That is, p = i^g (which is p= onpoints(i, g))
*    means that i is the index at slot p in the configuration g.
*    Equivalently, i = p^ig, where ig is the inverse of g.
*  - That convention is precisely the opposite to the one chosen by
*    Renato. Everything here and in xPerm follow Renato's convention.
*    There are two InversePerm actions in ToCanonicalOne in xTensor to
*    perform the change of notation. Here canonical_perm also has two
*    inverse  actions for backwards compatibility.
*  - The base of a SGS for the group S represents an ordering of the
*    slots of the tensor. Changing the base can be understood as a
*    change of priorities for index canonicalization.
*/

/**********************************************************************/

/* coset_rep. JMM, 30 June 2003
*
* This algorithm is encoded from Renato Portugal et al.
*
* This function gives a canonical representant of a n-permutation p
* in the group described by a SGS (pair base,GS) and the result is
* stored in cr. In other words, it gives the canonical representant
* of the right coset S.p. The free slots are different at return time.
*/

void coset_rep(int *p, int n,
               int *base, int bl, int *GS, int *m,
               int *freeps, int fl,
               int *cr)
	{

#ifdef VERBOSE_COSET						/*PPC*/
	printf("***** RIGHT-COSET-REP ALGORITHM ****\n");	/*PPC*/
	printf("Permutation ");
	print_perm(p, n, 1);		/*PPC*/
	printf("Base: ");
	print_list(base, bl, 1);		/*PPC*/
#endif								/*PPC*/
	/* Trivial case without symmetries */
	if (*m==0) {
		copy_list(p, cr, n);
		return;
		}
	/* else */

	int i, j, k, b, pp;
	int *deltap=       (int*)malloc(  n*sizeof(int)), deltapl;
	int *deltapsorted= (int*)malloc(  n*sizeof(int));
	int *om=           (int*)malloc(  n*sizeof(int));
	int *PERM=         (int*)malloc(  n*sizeof(int));
	int *perm2=        (int*)malloc(  n*sizeof(int));
	int *orbit=        (int*)malloc(  n*sizeof(int)), ol;
	int *orbit1=       (int*)malloc(  n*sizeof(int)), o1l;
	int *w=            (int*)malloc(  n*sizeof(int));
	int *nu=           (int*)malloc(n*n*sizeof(int));
	int *genset=       (int*)malloc(*m*n*sizeof(int)), gensetl;
	int *stab=         (int*)malloc(*m*n*sizeof(int)), mm;

	/* Copy p to PERM and GS to genset, to avoid side effects. */
	copy_list(p, PERM, n);
	copy_list(GS, genset, *m*n);
	gensetl=*m;

	/* Loop over elements of base */
	for(i=0; i<bl; i++) {
		b = base[i]; /* Slot to analyze */
		one_schreier_orbit(b, genset, gensetl, n,
		                   orbit, &ol, nu, w, 1);
#ifdef VERBOSE_COSET						/*PPC*/
		printf("\nAnalyzing slot %d\n", b);			/*PPC*/
		printf("orbit: ");
		print_list(orbit, ol, 1);		/*PPC*/
		printf("freeps: ");
		print_list(freeps, fl, 1);		/*PPC*/
#endif								/*PPC*/
		intersection(orbit, ol, freeps, fl, orbit1, &o1l);
#ifdef VERBOSE_COSET						/*PPC*/
		printf("Free slots that can go to that slot: ");	/*PPC*/
		print_list(orbit1, o1l, 1);				/*PPC*/
#endif								/*PPC*/
		if (o1l==0) continue; /* Slot with no symmetries */
		/* else */
		for(j=0; j<o1l; j++) {
			deltap[j] = onpoints(orbit1[j], PERM, n);
			}
		deltapl=o1l;
#ifdef VERBOSE_COSET						/*PPC*/
		printf("At those slots we resp. find indices deltap: ");/*PPC*/
		print_list(deltap, deltapl, 1);			/*PPC*/
#endif								/*PPC*/
		sortB(deltap, deltapsorted, deltapl, base, bl);
		k = position(deltapsorted[0], deltap, deltapl);
#ifdef VERBOSE_COSET						/*PPC*/
		printf("Least index: %d, at position k: %d of deltap\n",/*PPC*/
		       deltap[k-1], k);			/*PPC*/
#endif								/*PPC*/
		pp = orbit1[k-1];
#ifdef VERBOSE_COSET						/*PPC*/
		printf("That index is at tensor slot pp: %d\n", pp);	/*PPC*/
#endif								/*PPC*/
		/* Compute permutation om such that b^om = pp */
		trace_schreier(pp, nu, w, om, n);
#ifdef VERBOSE_COSET						/*PPC*/
		printf("We can move slot %d to slot %d with perm om:",	/*PPC*/
		       pp, b);					/*PPC*/
		print_perm(om, n, 0);
		printf(" in S\n");	/*PPC*/
#endif								/*PPC*/
		product(om, PERM, perm2, n);
#ifdef VERBOSE_COSET						/*PPC*/
		printf("New list of indices: ");			/*PPC*/
		print_perm(perm2, n, 1);			/*PPC*/
#endif								/*PPC*/
		copy_list(perm2, PERM, n);
		/* New positions of free indices */
		inverse(om, perm2, n);
		for(j=0; j<fl; j++) {
			freeps[j] = onpoints(freeps[j], perm2, n);
			}
#ifdef VERBOSE_COSET						/*PPC*/
		printf("Removing those perms that move slot %d\n", b);	/*PPC*/
#endif								/*PPC*/
		/* Note that we do not change base to have i as first
			member of base. This is not general, but I think
		here it is valid due to the nature of our problem */
		stabilizer(base+i, 1, genset, gensetl, n, stab, &mm);
		copy_list(stab, genset, mm*n);
		gensetl=mm;
		}
	/* Move to result */
	copy_list(PERM, cr, n);
	/* Return new SGS */
	copy_list(genset, GS, gensetl*n);
	*m=gensetl;
#ifdef VERBOSE_COSET						/*PPC*/
	printf("************ END OF ALGORITHM ***********\n");	/*PPC*/
#endif								/*PPC*/

	/* Free allocated memory */
	free(deltap);
	free(deltapsorted);
	free(om);
	free(PERM);
	free(perm2);
	free(orbit);
	free(orbit1);
	free(w);
	free(nu);
	free(genset);
	free(stab);

	}

/**********************************************************************/

/* SGSD. JMM, 26 June 2007
*
* We construct a number of functions that generate the strong
* generating set for a given list of (multiple) dummysets (all pairs
* of dummies of a vbundle) and repeatedsets (the positions of repeated
* indices, like components or directions). The main function is SGSD,
* but there are other elementary operations: SGSofdummyset,
* SGSofrepeatedset, movedummyset, moverepeatedset.
* We introduce another funny name here: drummy: either a dummy or a
* repeated index.
*/

/* SGS for a dummyset. List dummies not modified */
void SGSofdummyset(int *dummies, int dl, int sym, int n,
                   int *KD, int *KDl, int *bD, int *bDl)
	{

	if (dl==0) {
		*KDl=0;
		*bDl=0;
		return;
		} /* else */

	/* Number of pairs of dummies: dpl */
	int dpl = dl/2;
	int *range_perm = (int*)malloc(    n*sizeof(int));
	int *KD1 =        (int*)malloc(dpl*n*sizeof(int));
	int *KD2 =        (int*)malloc(dpl*n*sizeof(int));
	int i;

	range(range_perm, n);
	/* KD1: exchange indices. Ex: Cycles[{1,3},{2,4}]
	     There are always dpl-1 permutations */
	for (i=0; i<dpl-1; i++) {
		/* Copy the identity */
		copy_list(range_perm, KD1+i*n, n);
		/* Swap elements of consecutive pairs */
		KD1[ i*n+dummies[2*i  ]-1 ] = dummies[2*i+2];
		KD1[ i*n+dummies[2*i+2]-1 ] = dummies[2*i  ];
		KD1[ i*n+dummies[2*i+1]-1 ] = dummies[2*i+3];
		KD1[ i*n+dummies[2*i+3]-1 ] = dummies[2*i+1];
		}
	/* KD2: exchange indices of the same pair.
	     If sym!=0 there are dpl permutations */
	if (sym == 1) { /* Symmetric metric */
		for (i=0; i<dpl; i++) {
			/* Copy the identity */
			copy_list(range_perm, KD2+i*n, n);
			/* Swap elements of pair */
			KD2[i*n+dummies[2*i]-1] = dummies[2*i+1];
			KD2[i*n+dummies[2*i+1]-1] = dummies[2*i];
			/* Keep positive sign */
			}
		*KDl = 2*dpl-1;
		}
	else if (sym == -1) {
		for (i=0; i<dpl; i++) {
			/* Copy the identity */
			copy_list(range_perm, KD2+i*n, n);
			/* Swap elements of pair */
			KD2[i*n+dummies[2*i]-1] = dummies[2*i+1];
			KD2[i*n+dummies[2*i+1]-1] = dummies[2*i];
			/* Set negative sign */
			KD2[i*n+n-2] = n;
			KD2[i*n+n-1] = n-1;
			}
		*KDl = 2*dpl-1;
		}
	else if (sym == 0) {
		/* Do nothing */
		*KDl = dpl-1;
		}
	else {
		/* Unknown value of sym */
		}
	/* KD */
	copy_list(KD1, KD, (dpl-1)*n);
	if (sym!=0) copy_list(KD2, KD+(dpl-1)*n, dpl*n);
	/* base of group D */
	for (i=0; i<dpl; i++) {
		bD[i] = dummies[2*i];
		}
	*bDl=dpl;
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("KD: ");
	print_array_perm(KD, *KDl, n, 1);	/*PPC*/
	printf("bD: ");
	print_list(bD, *bDl, 1);		/*PPC*/
#endif								/*PPC*/
	free(range_perm);
	free(KD1);
	free(KD2);
	}

/* SGS for a repeatedset. List repes not modified */
void SGSofrepeatedset(int *repes, int rl, int n,
                      int *KD, int *KDl, int *bD, int *bDl)
	{

#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("From SGSofrepeatedset:\n");			/*PPC*/
	printf("repes: ");
	print_list(repes, rl, 1);
#endif								/*PPC*/
	if (rl==0) {
		*KDl=0;
		*bDl=0;
		return;
		} /* else */

	int *range_perm = (int*)malloc(    n*sizeof(int));
	int i;

	range(range_perm, n);
	for (i=0; i<rl-1; i++) {
		/* Copy the identity */
		copy_list(range_perm, KD+i*n, n);
		/* Swap elements of pair */
		KD[i*n+repes[i]-1] = repes[i+1];
		KD[i*n+repes[i+1]-1] = repes[i];
		/* Keep positive sign */
		}
	*KDl = rl-1;
	copy_list(repes, bD, rl-1);
	*bDl = rl-1;
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("KD: ");
	print_array_perm(KD, *KDl, n, 1);	/*PPC*/
	printf("bD: ");
	print_list(bD, *bDl, 1);		/*PPC*/
#endif								/*PPC*/
	free(range_perm);
	}

/* Move index in a dummyset. List dummies reordered */
void movedummyset(int firstd, int *dummies, int dl, int)
	{

	/* Find position of dummy and relative
		position of its pair */
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("Rearrange dummies for %d. dummies: ", firstd);	/*PPC*/
	print_list(dummies, dl, 1);				/*PPC*/
#endif								/*PPC*/
	int pos, j;
	pos = position(firstd, dummies, dl)-1;
	if (pos==-1) { /* firstd not in dummies */
		/* Do nothing */
		}
	else {
		int tmp;
		/* Swap all pairs if firstd found as second */
		if (pos%2==0) {
			/* 1st member: Do nothing */
			}
		else {
			pos=pos-1;
			/* 2nd member: Swap all pairs */
			for (j=0; j<dl/2; j++) {
				tmp= dummies[2*j];
				dummies[2*j]=dummies[2*j+1];
				dummies[2*j+1]=tmp;
				}
			}
		/* Exchange position of pair with first pair */
		if (pos==0) {
			/* 1st pair: Do nothing */
			}
		else {
			/* Exchange two pairs of dummies */
			tmp=dummies[0];
			dummies[0]=firstd;
			dummies[pos]=tmp;
			tmp=dummies[1];
			dummies[1]=dummies[pos+1];
			dummies[pos+1]=tmp;
			}
		}
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("Now dummies: ");                        	/*PPC*/
	print_list(dummies, dl, 1);				/*PPC*/
#endif								/*PPC*/
	}

/* Move index in a repeated set. List repes reordered */
void moverepeatedset(int firstd, int *repes, int rl)
	{

	/* Find position of dummy and relative
		position of its pair */
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("Rearrange dummies for %d. repes: ", firstd);	/*PPC*/
	print_list(repes, rl, 1);				/*PPC*/
#endif								/*PPC*/
	int pos;
	pos = position(firstd, repes, rl)-1;
	if (pos==-1) { /* firstd not in repes */
		/* Do nothing */
		}
	else {
		if (pos==0) {
			/* 1st index: Do nothing */
			}
		else {
			/* Exchange two positions */
			repes[pos] = repes[0];
			repes[0] = firstd;
			}
		}
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("Now repes: ");	                        	/*PPC*/
	print_list(repes, rl, 1);				/*PPC*/
#endif								/*PPC*/
	}

/* Remove a first-pair from dummies */
void dropdummyset(int firstd,
                  int *vds, int vdsl, int *dummies, int *dl)
	{

	int i, j, itotal=0;

	for (i=0; i<vdsl; i++) {
		if (dummies[itotal]==firstd && vds[i]!=0) {
			for (j=itotal; j<*dl-2; j++) {
				dummies[j] = dummies[j+2];
				}
			vds[i] = vds[i]-2;
			*dl = *dl - 2;
			return;
			}
		itotal = itotal + vds[i];
		}
	}

/* Remove a first-point from repes */
void droprepeatedset(int firstd,
                     int *vrs, int vrsl, int *repes, int *rl)
	{

	int i, j, itotal=0;

	for (i=0; i<vrsl; i++) {
		if (repes[itotal]==firstd && vrs[i]!=0) {
			for (j=itotal; j<*rl; j++) {
				repes[j] = repes[j+1];
				}
			vrs[i] = vrs[i]-1;
			*rl = *rl - 1;
			return;
			}
		itotal = itotal + vrs[i];
		}
	}

/* SGSD for a given list of dummies and repes */
void SGSD(int *vds, int vdsl, int *dummies, int dl, int *mQ,
          int *vrs, int vrsl, int *repes, int rl, int n,
          int firstd, int *KD, int *KDl, int *bD, int *bDl)
	{

	if (dl==0 && rl==0) {
		*KDl=0;
		*bDl=0;
		return;
		} /* else */

	int *tmp, tmpl;
	int *tmpGS   = (int*)malloc(n*n*sizeof(int)), tmpGSl;
	int *tmpbase = (int*)malloc(  n*sizeof(int)), tmpbasel;
	int i, itotal;

	/* Loop over all dummysets */
	itotal = 0;
	*KDl = 0;
	*bDl = 0;
	for (i=0; i<vdsl; i++) {
		tmpl = vds[i];
		tmp = dummies + itotal;
		movedummyset(firstd, tmp, tmpl, mQ[i]);
		itotal = itotal + tmpl;
		SGSofdummyset(tmp, tmpl, mQ[i], n,
		              tmpGS, &tmpGSl, tmpbase, &tmpbasel);
		copy_list(tmpGS, KD+(*KDl)*n, tmpGSl*n);
		*KDl = *KDl + tmpGSl;
		copy_list(tmpbase, bD+ *bDl, tmpbasel);
		*bDl = *bDl + tmpbasel;
		}

	/* Loop over all repeatedsets */
	itotal=0;
	for (i=0; i<vrsl; i++) {
		tmpl = vrs[i];
		tmp = repes + itotal;
		moverepeatedset(firstd, tmp, tmpl);
		itotal = itotal + tmpl;
		SGSofrepeatedset(tmp, tmpl, n,
		                 tmpGS, &tmpGSl, tmpbase, &tmpbasel);
		copy_list(tmpGS, KD+(*KDl)*n, tmpGSl*n);
		*KDl = *KDl + tmpGSl;
		copy_list(tmpbase, bD+ *bDl, tmpbasel);
		*bDl = *bDl + tmpbasel;
		}

	free(tmpGS);
	free(tmpbase);

#ifdef VERBOSE_DOUBLE
	printf("base of D:");
	print_list(bD, *bDl, 1);		/*PPC*/
	printf("GS of D (%d perms of degree %d):", *KDl, n);	/*PPC*/
	print_array_perm(KD, *KDl, n, 1);			/*PPC*/
#endif								/*PPC*/
	}

/**********************************************************************/

/* double_coset_rep. JMM, 1-5 July 2003.
* JMM, 9 November 2005. Corrected treatment of SGS of group D.
* JMM, 26 June 2007. Large extension to consider several dummysets and
*      several repeatedsets.
*
* This algorithm is encoded from Renato Portugal et al, with the
* extensions to repeatedsets.
*
* This function gives a canonical representant for the n-permutation g
* in the double coset S.g.D given by the groups S, described by a SGS
* (pair base, GS) and the group D, described by the direct product of
* the pair symmetric dpl pairs of dummies and the S_k groups of k
* repeated indices. Each dummyset has a metric_symmetry sign: if mQ=1
* the metric is symmetric. If mQ=-1 the metric is antisymmetric
* (spinor calculus). If mQ=0 there is no metric.
* Note that n = dl + dr. The result is stored in dcr.
*/

class alphastruct {
	public:
		alphastruct();
		alphastruct(const alphastruct&);
		~alphastruct();

		void init(int n);

		int *L; 	/* We assume that elements of L cannot
						be repeated. I only have experimental
							evidence for this */
		int Ll;
		int *s;
		int *d;
		int *o;

		int nn;
	};

alphastruct::alphastruct()
	: L(0), Ll(0), s(0), d(0), o(0), nn(0)
	{
	}

alphastruct::alphastruct(const alphastruct& other)
	{
	init(other.nn);
	for(int i=0; i<nn; ++i) {
		L[i]=other.L[i];
		s[i]=other.s[i];
		d[i]=other.d[i];
		o[i]=other.o[i];
		}
	Ll=other.Ll;
	}


void alphastruct::init(int n)
	{
	L = new int[n];
	s = new int[n];
	d = new int[n];
	o = new int[n];
	nn=n;
	}

alphastruct::~alphastruct()
	{
	if(L) delete [] L;
	if(s) delete [] s;
	if(d) delete [] d;
	if(o) delete [] o;
	nn=0;
	}


void TAB(std::vector<alphastruct>& ALPHA, int *L, int Ll, int *s1, int *d1, int n)
	{
	int i;
	int l=0;

	/* Search ALPHA for the l corresponding to L */
	//	std::cout << ALPHA.size() << ", " << Ll << std::endl;
	for (i=0; i<Ll; i++) l=ALPHA[l].o[L[i]-1];
	/* Copy permutations of element l */
	copy_list(ALPHA[l].s, s1, n);
	copy_list(ALPHA[l].d, d1, n);
	}

void F1(std::vector<alphastruct>& ALPHA, int *L, int Ll, int *g, int *list, int *listl, int n, int Deltabl, int *Deltab, int *DeltaD)
	{
	int c, c1, c2;
	int *sgd=  (int*)malloc(n*sizeof(int));
	int *TAB1= (int*)malloc(n*sizeof(int));
	int *TAB2= (int*)malloc(n*sizeof(int));
	int *tmp=  (int*)malloc(n*sizeof(int));

	TAB(ALPHA, L, Ll, TAB1, TAB2, n);

	/* Compute s.g.d */
	F2(TAB1, g, TAB2, sgd, n);
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("With L=");
	print_list(L, Ll, 0);		/*PPC*/
	printf(" we get sgd: ");
	print_perm(sgd, n, 1);		/*PPC*/
#endif								/*PPC*/
	/* Images of Deltab under sgd. Note that tmp has length
	Deltabl */
	for (c=0; c<Deltabl; c++)
		tmp[c] = onpoints(Deltab[c], sgd, n);
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf(" which maps slots in Deltab to indices list: ");/*PPC*/
	print_list(tmp, Deltabl, 1);			/*PPC*/
#endif								/*PPC*/
	/* Orbits of DeltaD containing the points in tmp */
	int oi;
	for (c1=0; c1<Deltabl; c1++) {
		oi = DeltaD[tmp[c1]-1];
		if(oi) {
			for(c2=0; c2<n; c2++) {
				if ((DeltaD[c2]==oi) &&
				      (!position(c2+1, list, *listl)))
					list[(*listl)++] = c2+1;
				}
			}
		}
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf(" whose points belong to orbits ");		/*PPC*/
	print_list(list, *listl, 1);				/*PPC*/
#endif								/*PPC*/
	free(sgd);
	free(TAB1);
	free(TAB2);
	free(tmp);
	}


/* Consistency check */
int consistency(int *array, int m, int n)
	{
	int *arrayp= (int*)malloc(m*n*sizeof(int)), arraypl;
	int *arrayn= (int*)malloc(m*n*sizeof(int)), arraynl;
	int i, ip, in, ret;

#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("Checking consistency in m:%d, n:%d\n", m, n);	/*PPC*/
	print_array_perm(array, m, n, 1);			/*PPC*/
#endif								/*PPC*/

	/* Detect sign of permutation */
	arraypl=0;
	arraynl=0;
	for(i=0; i<m; i++) {
		if (array[i*n+n-2]<array[i*n+n-1]) /* Positive */
			copy_list(array+i*n, arrayp+(arraypl++)*n, n);
		else                               /* Negative */
			copy_list(array+i*n, arrayn+(arraynl++)*n, n);
		}
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("Found positive perms: %d\n", arraypl);		/*PPC*/
	printf("Found negative perms: %d\n", arraynl);		/*PPC*/
#endif								/*PPC*/
	/* Here there are arraynl*arraypl comparisons. This
		should be improved with a better intersection
	algorithm which sorts the lists in advance */
	ret = 1; /* True */
	for (in=0; in<arraynl; in++) {
		for (ip=0; ip<arraypl; ip++) {
			if (equal_list(arrayp+ip*n, arrayn+in*n, n-2)) {
				ret = 0; /* False */
				break;
				}
			}
		}
#ifdef VERBOSE_DOUBLE						/*PPC*/
	if (ret) printf("Found no problem in check\n");		/*PPC*/
	else printf("Found perm with two signs.\n");		/*PPC*/
#endif								/*PPC*/
	free(arrayp);
	free(arrayn);
	return(ret);
	}


void double_coset_rep(int *g, int n, int *base, int bl, int *GS, int m,
                      int *vds, int vdsl, int *dummies, int dl, int *mQ,
                      int *vrs, int vrsl, int *repes, int rl, int *dcr)
	{

	int i, j, l, jj, kk, c;         /* Counters */
	int result;
	/* Inverse of permutation g */
	int *ig=            (int*)malloc(  n*sizeof(int));
	/* All drummies, both the pair-dummies and the repes */
	int *drummies=      (int*)malloc(  n*sizeof(int)), dril;
	/* The initial slots of all those dummies */
	int *drummyslots=   (int*)malloc(  n*sizeof(int));
	/* Temporary space for sorting */
	int *drummytmp=     (int*)malloc(  n*sizeof(int)), drummytmpl;
	int *drummytmp2=    (int*)malloc(  n*sizeof(int));
	/* Bases for group S */
	int *bS=            (int*)malloc(  n*sizeof(int)), bSl;
	int *bSsort=        (int*)malloc(  n*sizeof(int));
	/* gs for group S. It cannot be larger than GS */
	int *KS=            (int*)malloc(m*n*sizeof(int)), KSl;
	/* gs for group D. The number dl+dr is an upper bound */
	int *KD=            (int*)malloc(n*n*sizeof(int)), KDl;
	int *bD=            (int*)malloc(  n*sizeof(int)), bDl;
	int *nu=            (int*)malloc(n*n*sizeof(int));
	int *w=             (int*)malloc(  n*sizeof(int));
	int *Deltab=        (int*)malloc(  n*sizeof(int)), Deltabl;
	int *DeltaD=        (int*)malloc(  n*sizeof(int));
	int *IMAGES=        (int*)malloc(  n*sizeof(int)), IMAGESl;
	int *IMAGESsorted=  (int*)malloc(  n*sizeof(int));
	int *p=             (int*)malloc(  n*sizeof(int));
	int *nuD=           (int*)malloc(n*n*sizeof(int));
	int *wD=            (int*)malloc(  n*sizeof(int));
	int *Deltap=        (int*)malloc(  n*sizeof(int)), Deltapl;
	int *NEXT=          (int*)malloc(  n*sizeof(int)), NEXTl;
	int *L=             (int*)malloc(  n*sizeof(int)), Ll;
	int *L1=            (int*)malloc(  n*sizeof(int)), L1l;
	int *s=             (int*)malloc(  n*sizeof(int));
	int *d=             (int*)malloc(  n*sizeof(int));
	int *list1=         (int*)malloc(  n*sizeof(int)), list1l;
	int *list2=         (int*)malloc(  n*sizeof(int)), list2l;
	int *perm1=         (int*)malloc(  n*sizeof(int));
	int *perm2=         (int*)malloc(  n*sizeof(int));
	int *perm3=         (int*)malloc(  n*sizeof(int));
	int *s1=            (int*)malloc(  n*sizeof(int));
	int *d1=            (int*)malloc(  n*sizeof(int));


	/* We use Renato's notation, with g mapping slots to indices */

	/**************************************************************
	*                       DEFINITIONS
	**************************************************************/

	/* Given a list L={i1, ..., il} the function TAB must return
	a pair (s1, d1) corresponding to L. The whole collection of
		lists L and their pairs is stored in the array ALPHA of
		structures alphastruct. Each L in ALPHA is actually
	identified by its position l in the array. The structure
	for a given L contains a list of values of l corresponding
	to the n possible extensions of L.
	*/

	/* Define ALPHA and TAB */

	int ALPHAl;
	int *ALPHAstep= (int*)malloc(n*sizeof(int));

	/* Initialize ALPHA to {} and TAB to {id, id} */
	ALPHAl= 1;
	std::vector<alphastruct> ALPHA(1);
	//	alphastruct *ALPHA= new alphastruct[1]; //(struct alphastruct*)malloc(sizeof(struct alphastruct));
	ALPHA[0].init(n);
	ALPHA[0].Ll=0;
	Ll=0; // ?
	range(ALPHA[0].s, n);
	range(ALPHA[0].d, n);
	ALPHAstep[0]=0;
	ALPHAstep[1]=1;


	/**************************************************************
	*               CONSTRUCTION OF BASES
	**************************************************************/

	/* Join all drummies */
	copy_list(dummies, drummies, dl);
	copy_list(repes, drummies+dl, rl);
	dril = dl + rl;

	/* Slots of all those drummies */
	inverse(g, ig, n);
	for (i=0; i<dril; i++) {
		drummyslots[i] = onpoints(drummies[i], ig, n);
		}

	/* Initialize KS */
	copy_list(GS, KS, m*n);
	KSl=m;

	/* Extend base to bS to cover all positions of drummies.
	* We assume that in the intersection we get bS=base really.
	* We sort the missing drummies, but not the given base  */
	intersection(base, bl, drummyslots, dril, bS, &bSl);
	complement(drummyslots,dril, base,bl, 1, drummytmp,&drummytmpl);
	sort(drummytmp, drummytmp2, drummytmpl);
	copy_list(drummytmp2, bS+bSl, drummytmpl);
	bSl = bSl + drummytmpl;
#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("All drummies: drummies: ");			/*PPC*/
	print_list(drummies, dril, 1);				/*PPC*/
	printf("Their slots: drummyslots: ");			/*PPC*/
	print_list(drummyslots, dril, 1);			/*PPC*/
	printf("base: ");					/*PPC*/
	print_list(base, bl, 1);				/*PPC*/
	printf("Complement: drummytmp: ");			/*PPC*/
	print_list(drummytmp, drummytmpl, 1);			/*PPC*/
	printf("Sort them: drummytmp2: ");			/*PPC*/
	print_list(drummytmp2, drummytmpl, 1);			/*PPC*/
	printf("base extended to bS: ");			/*PPC*/
	print_list(bS, bSl, 1);					/*PPC*/
#endif								/*PPC*/
	/* Generate associated base for sorting names of dummies */
	/* We choose a particular form of bSsort for aesthetics */
	sort(drummies,drummytmp2,dril);
	sort(bS, drummytmp, bSl);
	for (i=0; i<bSl; i++) {
		bSsort[i]=drummytmp2[position(bS[i],drummytmp,bSl)-1];
		}

#ifdef VERBOSE_DOUBLE						/*PPC*/
	printf("drummies: ");					/*PPC*/
	print_list(drummies, dril, 1);				/*PPC*/
	printf("Base of SGSS: bS: ");				/*PPC*/
	print_list(bS, bSl, 1);					/*PPC*/
	printf("Sorted slots: ");				/*PPC*/
	print_list(drummytmp2, bSl, 1);				/*PPC*/
	printf("Sorted base bS: ");				/*PPC*/
	print_list(drummytmp, bSl, 1);				/*PPC*/
	printf("Base for sorting: bSsort: ");			/*PPC*/
	print_list(bSsort, bSl, 1);				/*PPC*/
#endif								/*PPC*/


	/*******************/
	/**** Main loop ****/
	/*******************/

	/* Note we use i=1..bSl instead of the usual i=0 ..(bSl-1) */
	for (i=1; i<=bSl; i++) {
		int b=bS[i-1];

#ifdef VERBOSE_DOUBLE						/*PPC*/
		printf("\n************** Loop i=%d *************\n", i);/*PPC*/
		printf("Analyzing slot bS[%d]=%d of tensor\n", i-1, b);	/*PPC*/
#endif								/*PPC*/
		/* Schreier vector of S */
		schreier_vector(b, KS, KSl, n, nu, w);
		one_orbit(b, KS, KSl, n, Deltab, &Deltabl);
#ifdef VERBOSE_DOUBLE						/*PPC*/
		printf("Under S, slot %d go to slots Deltab: ",b);	/*PPC*/
		print_list(Deltab, Deltabl, 1);				/*PPC*/
#endif								/*PPC*/
		/* Compute SGS for group D. Do not rearrange drummies */
		SGSD(vds, vdsl, dummies, dl, mQ,
		     vrs, vrsl, repes, rl, n,
		     0, KD, &KDl, bD, &bDl);
		/* Orbits of D */
		all_orbits(KD, KDl, n, DeltaD);
#ifdef VERBOSE_DOUBLE						/*PPC*/
		printf("Orbits of indices: DeltaD: ");			/*PPC*/
		print_list(DeltaD, n, 1);				/*PPC*/
#endif								/*PPC*/
		/* Images of b under elements of S.g.D.
		Deltab and DeltaD are used by F1 */
		IMAGESl=0;
		for (c=ALPHAstep[i-1]; c<ALPHAstep[i]; c++)
			F1(ALPHA, ALPHA[c].L, ALPHA[c].Ll, g, IMAGES, &IMAGESl, n, Deltabl, Deltab, DeltaD);
#ifdef VERBOSE_DOUBLE						/*PPC*/
		printf("At slot %d we can have indices IMAGES: ", b);	/*PPC*/
		print_list(IMAGES, IMAGESl, 1);				/*PPC*/
		printf("IMAGESl: %d\n", IMAGESl);			/*PPC*/
#endif								/*PPC*/
		/* If there are no images we have finished */
		if(IMAGESl==0) continue;

		/* Find minimum index */
		sortB(IMAGES, IMAGESsorted, IMAGESl, bSsort, bSl);
		p[i-1] = IMAGESsorted[0];
#ifdef VERBOSE_DOUBLE						/*PPC*/
		printf("The least of them is p[i-1]=%d\n", p[i-1]);	/*PPC*/
#endif								/*PPC*/
		/* Recompute SGS of D */
		if (dl>0 || rl>0) {
			SGSD(vds, vdsl, dummies, dl, mQ,
			     vrs, vrsl, repes, rl, n,
			     p[i-1], KD, &KDl, bD, &bDl);
			}
		else {   /* Do nothing */
#ifdef VERBOSE_DOUBLE						/*PPC*/
			printf("Rearrangement of base of D not required.\n");	/*PPC*/
#endif								/*PPC*/
			}
		/* Schreier vector of D */
		schreier_vector(p[i-1], KD, KDl, n, nuD, wD);
		/* Orbit of p[i-1] under D */
		one_orbit(p[i-1], KD, KDl, n, Deltap, &Deltapl);
#ifdef VERBOSE_DOUBLE						/*PPC*/
		printf("The orbit of index %d is Deltap: ", p[i-1]);	/*PPC*/
		print_list(Deltap, Deltapl, 1);				/*PPC*/
		printf("Looking for digs moving index %d to slot %d\n",	/*PPC*/
		       p[i-1], b);					/*PPC*/
#endif								/*PPC*/
		/* Calculate ALPHA and TAB */
		ALPHAstep[i+1]=ALPHAstep[i];
		for(l=ALPHAstep[i-1]; l<ALPHAstep[i]; l++) {
#ifdef VERBOSE_DOUBLE						/*PPC*/
			printf("Loop with l=%d\n", l);				/*PPC*/
#endif								/*PPC*/
			Ll=ALPHA[l].Ll;
			copy_list(ALPHA[l].L, L, Ll);
#ifdef VERBOSE_DOUBLE						/*PPC*/
			printf("L: ");
			print_list(L, Ll, 1);			/*PPC*/
#endif								/*PPC*/
			copy_list(ALPHA[l].s, s, n);
			copy_list(ALPHA[l].d, d, n);
#ifdef VERBOSE_DOUBLE						/*PPC*/
			printf("TAB[L]={");					/*PPC*/
			print_perm(s, n, 0);					/*PPC*/
			printf(", ");						/*PPC*/
			print_perm(d, n, 0);					/*PPC*/
			printf("}\n");						/*PPC*/
#endif								/*PPC*/
			list1l=Deltabl;
			for (c=0; c<list1l; c++)
				list1[c]=onpoints(Deltab[c], s, n);
			product(g, d, perm1, n);
			inverse(perm1, perm2, n);
			list2l=Deltapl;
			for (c=0; c<list2l; c++)
				list2[c]=onpoints(Deltap[c], perm2, n);
#ifdef VERBOSE_DOUBLE						/*PPC*/
			printf("NEXT: intersection of sets of slots ");		/*PPC*/
			print_list(list1, list1l, 0);
			printf(" and ");		/*PPC*/
			print_list(list2, list2l, 1);				/*PPC*/
#endif								/*PPC*/
			intersection(list1, list1l, list2, list2l,
			             NEXT, &NEXTl);
#ifdef VERBOSE_DOUBLE						/*PPC*/
			printf("Intermediate slots NEXT: ");			/*PPC*/
			print_list(NEXT, NEXTl, 1);				/*PPC*/
#endif								/*PPC*/

			for(jj=0; jj<NEXTl; jj++) {
				j = NEXT[jj];
				inverse(s, perm1, n);
				trace_schreier(onpoints(j,perm1,n), nu,w, perm2,n);
				product(perm2, s, s1, n);
#ifdef VERBOSE_DOUBLE						/*PPC*/
				printf("From slot %d to interm. slot %d use s1: ",	/*PPC*/
				       b, j);
				print_perm(s1, n, 1);			/*PPC*/
#endif								/*PPC*/
				product(g, d, perm2, n);
				trace_schreier(onpoints(j,perm2,n),nuD,wD,perm3,n);
				inverse(perm3, perm1, n);
				product(d, perm1, d1, n);
				copy_list(L, L1, Ll);
				L1l=Ll;
				L1[L1l++] = j;
#ifdef VERBOSE_DOUBLE						/*PPC*/
				printf("d1: ");
				print_perm(d1, n, 1);			/*PPC*/
				printf("L1: ");
				print_list(L1, L1l, 1);			/*PPC*/
#endif								/*PPC*/

				kk=ALPHAstep[i+1];
				ALPHA[l].o[j-1] = kk;
				ALPHAl++;
				ALPHAstep[i+1]++;
				//				std::cout << "resizing " << ALPHAl << std::endl;
				//				ALPHA = (struct alphastruct*)realloc(ALPHA, ALPHAl* sizeof(struct alphastruct));
				ALPHA.resize(ALPHAl);
				ALPHA.back().init(n);
				copy_list(L1, ALPHA[kk].L, L1l);
				ALPHA[kk].Ll = L1l;
				copy_list(s1, ALPHA[kk].s, n);
				copy_list(d1, ALPHA[kk].d, n);

				F2(s1, g, d1, perm1, n);
#ifdef VERBOSE_DOUBLE						/*PPC*/
				printf("This gives the new index configuration: ");	/*PPC*/
				inverse(perm1, perm2, n);				/*PPC*/
				print_perm(perm2, n, 1);				/*PPC*/
				for(ii=0; ii<i; ii++) {					/*PPC*/
					product(s1, g, perm2, n);			/*PPC*/
					product(perm2, d1, perm3, n);			/*PPC*/
					printf("Checking slot %d with point %d: ",	/*PPC*/
					       bS[ii], p[ii]);				/*PPC*/
					if(onpoints(bS[ii], perm3, n)==p[ii])		/*PPC*/
						printf("True\n");			/*PPC*/
					else printf("*************FALSE************\n");/*PPC*/
					}							/*PPC*/
#endif								/*PPC*/
				}
			}
			{
			/* Verify if there are 2 equal permutations
			of opposite sign in SgD */
			int *array= (int*)malloc(n*(ALPHAstep[i+1]-ALPHAstep[i])
			                         *sizeof(int));
			int arrayl= 0;
#ifdef VERBOSE_DOUBLE						/*PPC*/
			printf("Astep[i-1]=%d, Astep[i]=%d, Astep[i+1]=%d\n",	/*PPC*/
			       ALPHAstep[i-1], ALPHAstep[i], ALPHAstep[i+1]);	/*PPC*/
#endif 								/*PPC*/
			for(l=ALPHAstep[i]; l<ALPHAstep[i+1]; l++) {
				F2(ALPHA[l].s, g, ALPHA[l].d, perm1, n);
				copy_list(perm1, array+(arrayl++)*n, n);
				}
#ifdef VERBOSE_DOUBLE						/*PPC*/
			printf("Perform check.\n");				/*PPC*/
#endif								/*PPC*/
			result= consistency(array, arrayl, n);
#ifdef VERBOSE_DOUBLE						/*PPC*/
			printf("Result of check: %d\n", result);		/*PPC*/
#endif								/*PPC*/
			free(array);
			if (!result) break;
			}
		/* Find the stabilizers S^(i+1) and D^(i+1) */
		stabilizer(&bS[i-1], 1, KS, KSl, n, KS, &KSl);
		dropdummyset(p[i-1], vds, vdsl, dummies, &dl);
		droprepeatedset(p[i-1], vrs, vrsl, repes, &rl);
#ifdef VERBOSE_DOUBLE						/*PPC*/
		printf("Remove perms of KS moving slot %d\n", b);	/*PPC*/
		printf("Remove perms of KD moving index %d\n", p[i-1]);	/*PPC*/
#endif								/*PPC*/

		} /* End of main loop */

	/* Result */
	if (result==0) {
		zeros(perm1, n);
		}
	else {
		l=ALPHAstep[i-1];
		F2(ALPHA[l].s, g, ALPHA[l].d, perm1, n);
		}
	copy_list(perm1, dcr, n);

	/* Free allocated memory */
	//	free(ALPHA);
	free(ALPHAstep);
	free(ig);
	free(drummies);
	free(drummyslots);
	free(drummytmp);
	free(drummytmp2);
	free(bS);
	free(bSsort);
	free(KS);
	free(KD);
	free(bD);
	free(nu);
	free(w);
	free(Deltab);
	free(DeltaD);
	free(IMAGES);
	free(IMAGESsorted);
	free(p);
	free(nuD);
	free(wD);
	free(Deltap);
	free(NEXT);
	free(L);
	free(L1);
	free(s);
	free(d);
	free(list1);
	free(list2);
	free(perm1);
	free(perm2);
	free(perm3);
	free(s1);
	free(d1);
	}

/**********************************************************************/

/* canonical_perm. JMM, 5 July 2003
*
* JMM, 26 June 2007. Function extended. The old function is kept
*      as a call to the new function for backwards compatibility.
*
* The "dummy" group D is given through the list dummyps of dpl pairs
* of (initial) slots of dummies. The list freeps (length fl) contains
* the slots of the free indices. Clearly 2*dpl+fl+2=n.
* See notes for canonical_perm_ext.
* Parameter ob is now useless. Kept for backwards compatibility.
*/

void canonical_perm(int *PERM,
                    int SGSQ, int *base, int bl, int *GS, int m, int n,
                    int *freeps, int fl, int *dummyps, int dpl, int, int metricQ,
                    int *CPERM)
	{

	int i;
	int vds;
	int mQ;
	int *repes= NULL;
	int *vrs= NULL;
	int *PERM1=   (int*)malloc(n    *sizeof(int));
	int *PERM2=   (int*)malloc(n    *sizeof(int));
	int *frees=   (int*)malloc(fl   *sizeof(int));
	int *dummies= (int*)malloc(2*dpl*sizeof(int));

	/* Construct "vectors" vds and mQ */
	vds = 2*dpl;
	mQ = metricQ;

	/* !!!!!!!! Change to Renato's notation !!!!!!!! */
	inverse(PERM, PERM1, n);
	for (i=0; i<fl; i++) {
		frees[i] = onpoints(freeps[i], PERM1, n);
		}
	for (i=0; i<2*dpl; i++) {
		dummies[i] = onpoints(dummyps[i], PERM1, n);
		}

	/* Call new, extended function */
	canonical_perm_ext(PERM1, n, SGSQ, base, bl, GS, m,
	                   frees, fl, &vds, 1, dummies, 2*dpl, &mQ,
	                   vrs, 0, repes, 0,
	                   PERM2);

	/* !!!!!!!! Change back to our notation !!!!!!!! */
	if (PERM2[0] != 0) inverse(PERM2, CPERM, n);
	else copy_list(PERM2, CPERM, n);

	/* Free allocated space */
	free(PERM1);
	free(PERM2);
	free(frees);
	free(dummies);
	}

/**********************************************************************/

/* canonical_perm_ext. JMM, 26 June 2007
*
* This function finds a canonical representant of the permutation PERM
* in the double_coset with group S (slot symmetries) and group D
* (dummy index or repeated index symmetries).
*
* There are two steps: first S, then S and D. The "slot" group S is
* given through the generating set GS (with m n-permutations), that
* can be a SGS (then SGSQ must be 1) or not (then SGSQ must be 0).
* In the former case base (length bl) contains the associated base;
* in the latter case a SGS will be constructed using the Schreier_Sims
* algorithm using the points in base as the first points for the base.
* The algorithm first calls coset_rep with PERM and the free indices
* converting PERM into PERM1. The SGS is then stabilized with respect
* to those points. Then the double_coset_rep algorithm is called with
* PERM1 returning PERM2, which is finally copied to CPERM.
*
* The "drummy" group D is given through the specification of so-called
* dummysets and repeatedsets. A dummyset is a collection of pairs of
* indices (first the up-index, then the down-index) and a symmetry
* switch saying if there is a symmetric metric (switch 1), an
* antisymmetric metric (-1) or no metric at all (0). A repeated set
* is simply a list of the names of repeated indices. Both dummies and
* repes are collectively called drummies. Note that dummysets and
* repeatedsets contain **names** (i.e. positions in the canonical
* configuration) in the initial configuration.
*
* The result is stored in the permutation CPERM.
* Note that we input a pointer to GS, and not a pointer to a pointer
* because there is no need to return a modified GS.
*
* This is probably the most important function of this code, and hence
* it is worth to explain the meaning of the input fields:
*      PERM: pointer to permutation (Images notation) to canonicalize
*      n:    degree of perm (length of list of images)
*      SGSQ: switch. 1 if we supply a SGS, 0 if we only give a GS
*      base: sorted list of points for the SGS. Start of base if SGSQ=0
*      bl:   length of base
*      GS:   pointer to the list of permutations forming the GS
*      m:    number of permutations in the GS
*      freeps: list of initial names of free indices
*      fl:   length of list freeps
*      vds:  list of lengths of dummysets
*      vdsl: length of vds (number of dummysets)
*      dummies: list with pairs of names dummies
*      dl:   length of list dummies (sum of elements of vdsl)
*      mQ:   list (of length vdsl) with symmetry signs
*      vrs:  list of lengths of repeatedsets
*      vrsl: length of rds (number of repeatedsets)
*      repes: list with repeated indices
*      rl:   length of list repes (sum of elements of vrsl)
*      CPERM: pointer to return the canonicalized permutation
*/

void canonical_perm_ext(int *PERM, int n,
                        int SGSQ, int *base, int bl, int *GS, int m,
                        int *frees, int fl,
                        int *vds, int vdsl, int *dummies, int dl, int *mQ,
                        int *vrs, int vrsl, int *repes, int rl,
                        int *CPERM)
	{

	int i;
	int *freeps=    (int*)malloc(fl*sizeof(int));
	int *PERM1=     (int*)malloc(n *sizeof(int));
	int *PERM2=     (int*)malloc(n *sizeof(int));
	int *newbase=   (int*)malloc(n *sizeof(int)), newbl;
	int **newGS = NULL;
	int *pointer;
	int newm;
	int *tmpbase=   (int*)malloc(n *sizeof(int)), tmpbl;
	int num=0;

	pointer= (int*)malloc(m*n*sizeof(int));
	newGS= &pointer;

#ifdef VERBOSE_CANON						/*PPC*/
	printf("can: input base: ");				/*PPC*/
	print_list(base, bl, 1);				/*PPC*/
#endif								/*PPC*/

	if (!SGSQ) { /* Compute a Strong Generating Set */
		nonstable_points(base, bl, GS, m, n,
		                 tmpbase, &tmpbl);
		schreier_sims(tmpbase, tmpbl, GS, m, n,
		              newbase, &newbl, newGS, &newm, &num);
		}
	else {
		copy_list(base, newbase, bl);
		newbl=bl;
		copy_list(GS, *newGS, m*n);
		newm=m;
		}

#ifdef VERBOSE_CANON						/*PPC*/
	printf("can: SGS computed.\n");				/*PPC*/
	printf("can: newbase: ");				/*PPC*/
	print_list(newbase, newbl, 1);				/*PPC*/
#endif								/*PPC*/

	/* Compute slots of free indices */
	inverse(PERM, PERM1, n);
	for (i=0; i<fl; i++) {
		freeps[i] = onpoints(frees[i], PERM1, n);
		}
	/* Call coset_rep algorithm. Result in PERM1 */
	coset_rep(PERM, n, newbase, newbl, *newGS, &newm,
	          freeps, fl, PERM1);
#ifdef VERBOSE_CANON						/*PPC*/
	printf("can: Canonical perm after coset algorithm: ");	/*PPC*/
	print_perm(PERM1, n, 1);				/*PPC*/
	printf("New positions of free indices: ");		/*PPC*/
	print_list(freeps, fl, 1);				/*PPC*/
#endif								/*PPC*/

	if (dl+rl==0) { /* No drummy indices */
		copy_list(PERM1, CPERM, n);
		}
	else {

		complement(newbase, newbl, freeps, fl, 1, tmpbase, &tmpbl);
		copy_list(tmpbase, newbase, tmpbl);
		newbl=tmpbl;
		stabilizer(freeps, fl, *newGS, newm, n, *newGS, &newm);
#ifdef VERBOSE_CANON						/*PPC*/
		printf("can: newbase after fixing: ");			/*PPC*/
		print_list(newbase, newbl, 1);				/*PPC*/
		printf("can: newGS after fixing: ");			/*PPC*/
		print_array_perm(*newGS, newm, n, 1);			/*PPC*/
#endif								/*PPC*/

		/* Apply dummy-indices algorithm. Result in PERM2 */
#ifdef VERBOSE_CANON						/*PPC*/
		printf("can: Starting double_coset algorithm.\n");	/*PPC*/
#endif								/*PPC*/
		double_coset_rep(PERM1, n, newbase, newbl, *newGS, newm,
		                 vds, vdsl, dummies, dl, mQ,
		                 vrs, vrsl, repes, rl, PERM2);
#ifdef VERBOSE_CANON						/*PPC*/
		printf("can: Finished double_coset algorithm.\n");	/*PPC*/
#endif								/*PPC*/

		/* Copy to result */
		copy_list(PERM2, CPERM, n);

		}

	/* Free allocated memory */
	free(freeps);
	free(PERM1);
	free(PERM2);
	free(newbase);
	free(tmpbase);
	free(*newGS);

#ifdef VERBOSE_CANON						/*PPC*/
	printf("************ END OF ALGORITHM ***********\n");	/*PPC*/
#endif								/*PPC*/
	}