/*********************************************************************
 *********************************************************************
 *********************************************************************
 *
 *  xperm.c
 *
 * 	(C) Jose M. Martin-Garcia 2003-2011.
 *      jose@xAct.es
 *      http://www.xAct.es/
 *
 *	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.
 *
 *  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
 *
 *********************************************************************
 *********************************************************************
 *********************************************************************/

/*********************************************************************
 *                             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);

void stab_chain(int *base, int bl, int *GS, int m, int n, int **chain, int *cl);

void one_orbit_chain(int point, int *GS, int *poslist, int poslistl, int n, int *orbit, int *ol);

void one_schreier_orbit_chain(int point, int *GS, int *poslist, int poslistl,
	int n, int *orbit, int *ol, int *nu, int *w, int init);

void conjugate_chain(int *base, int bl, int *GS, int m, int n, int *p);

void basechange_chain(int **base, int *bl, int **GS, int *m, int n, int ***chain, int **cl, int *newbase, int newbl);

void appendbasepoint_chain(int **base, int *bl, int ***chain, int **cl, int newbasepoint);

void interchange_chain(int **base, int *bl, int **GS, int *m, int n, int ***chain, int **cl, int j);


/*********************************************************************
 *                         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");
}

void fprint_perm(FILE *file, int *p, int n, int nl) {
        int i;
        if (isid(p,n)) fprintf(file, "id");
        else {
                fprintf(file,"(");
                fprintf(file,"%d", p[0]);
                for (i=1; i<n; i++) {
                        fprintf(file,",%d", p[i]);
                }
                fprintf(file,")");
        }
        if (nl) fprintf(file,"\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");
}

void fprint_list(FILE *file, int *list, int n, int nl) {
        int i;
        fprintf(file, "{");
        if (n>0) fprintf(file, "%d", list[0]);            /* No comma */
        for (i=1; i<n; i++) fprintf(file, ",%d", list[i]);/* Comma separated */
        fprintf(file,"}");
        if (nl) fprintf(file,"\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).
 */

}

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

/* 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 or=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] = or;
			/* Increment index for next orbit */
			or++;
		}
	}

	/* 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, 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 */
#ifdef VERBOSE_COSET						/*PPC*/
	printf("GS before stabilizing:\n");			/*PPC*/
	print_array_perm(genset, gensetl, n, 1);		/*PPC*/
#endif								/*PPC*/
		stabilizer(base+i, 1, genset, gensetl, n, stab, &mm);
#ifdef VERBOSE_COSET						/*PPC*/
	printf("GS after stabilizing:\n");			/*PPC*/
	print_array_perm(stab, mm, n, 1);			/*PPC*/
#endif								/*PPC*/
#ifdef VERBOSE_COSET						/*PPC*/
	printf("Removing those perms that move slot %d\n", b);	/*PPC*/
#endif								/*PPC*/
		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,j;

	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 sym) {

	/* 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.
 */

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, ii, jj, kk, c;         /* Counters */
	int oi, 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 */
	struct alphastruct {
		/* L */
		int L[n]; 	/* We assume that elements of L cannot
				   be repeated. I only have experimental
				   evidence for this */
		int Ll;
		/* s, d */
		int s[n];
		int d[n];
		/* other */
		int o[n];
				   
	} *ALPHA;
	int ALPHAl;
	int *ALPHAstep= (int*)malloc(n*sizeof(int));

	/* Initialize ALPHA to {} and TAB to {id, id} */
	ALPHAl= 1;
	ALPHA= (struct alphastruct*)malloc(sizeof(struct alphastruct));
	ALPHA[0].Ll=0;
	range(ALPHA[0].s, n);
	range(ALPHA[0].d, n);
	ALPHAstep[0]=0;
	ALPHAstep[1]=1;

	/* TAB */
	void TAB(int *L, int Ll, int *s1, int *d1, int n) {

		int i;
		int l=0;
		
		/* Search ALPHA for the l corresponding to L */
		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);

	} /* End of function TAB */

	/* Subroutines F1 and F2 */
	/* 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 */

	void F1(int *L, int Ll, int *g, int *list, int *listl, int n) {

		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(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 */
		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);
	} /* End of function F1 */

        /* 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);

	} /* End of funcion consistency */

        /**************************************************************
         *               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[c].L, ALPHA[c].Ll, g, IMAGES, &IMAGESl, n);
#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]++;
		    ALPHA = (struct alphastruct*)realloc(ALPHA, ALPHAl*
				    sizeof(struct alphastruct));
		    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 ob, 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*/
        printf("can: newGS:\n");				/*PPC*/
	print_array_perm(*newGS, m, n, 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;
#ifdef VERBOSE_CANON
	printf("can: newGS before fixing: ");			/*PPC*/
	print_array_perm(*newGS, newm, n, 1);			/*PPC*/
#endif
	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*/
}

/**********************************************************************
 *
 * Backtrack search
 *
 * We give a SGS and a property identifying a subgroup K of it.
 * We return a SGS for K with the same base.
 * So far properties can be:
 *      1. Centralizer of a permutation
 *      2. Normalizer of a subgroup
 *      3. Intersection of two subgroups
 *      4. SetStabilizer of a set
 */

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

        int *tmp  = (int *)malloc(n*sizeof(int)), tmpl;
        int *tmp2 = (int *)malloc(n*sizeof(int));

        one_orbit(point, GS, m, n, tmp, &tmpl);
        sortB(tmp, tmp2, tmpl, base, bl);
        one_schreier_orbit(tmp2[0],GS,m,n,orbit,ol,nu,w,init);

        free(tmp);
        free(tmp2);

}

void sorted_schreier_vector(int point, int *base, int bl, 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 */
        sorted_schreier_orbit(point, base, bl, GS, m, n, orbit, &ol, nu, w, 1);
        while(ol--) usedpoints[j++] = orbit[ol];
        /* Other orbits. Do not initialize vector */
        for(i=1; i<=n; i++) {
                if (!position(i, usedpoints, j)) {
                        sorted_schreier_orbit(i, base, bl, GS, m, n, orbit, &ol, nu, w, 0);
                        while(ol--) usedpoints[j++] = orbit[ol];
                }
        }
        /* Free allocated memory */
        free(orbit);
        free(usedpoints);

}


/* This function constructs the permutation corresponding to
   some list of images (with length ll<=bl). At the end the
   list is identical to the first ll points of base.
 */

void from_base_image(int *images, int ll,
        int *base, int bl, int *GS, int m, int n, int *perm) {

        int *list= (int *)malloc(bl*    sizeof(int));
        int *nu=   (int *)malloc(   n*n*sizeof(int));
        int *w=    (int *)malloc(     n*sizeof(int));
        int *u=    (int *)malloc(     n*sizeof(int));
        int *tmp=  (int *)malloc(     n*sizeof(int));
        int *stab= (int *)malloc(   m*n*sizeof(int)), sl;
        int i,j;

        copy_list(images,list,ll);
        range(perm,n);
        for (i=0; i<ll; i++) {
                /* It would be better to stabilize stabilizers,
                   instead of stabilizing the original set GS.
                   I think we are computing this stabilizers too
                   many times */
                stabilizer(base,i,GS,m,n,stab,&sl);
                sorted_schreier_vector(list[i], base+i, bl-i, stab, sl, n, nu, w);
                trace_schreier(list[i], nu, w, u, n);
                product(u,perm,tmp,n);
                copy_list(tmp,perm,n);
                inverse(u,tmp,n);
                for (j=0; j<ll; j++) {
                        list[j]=onpoints(list[j],tmp,n);
                }
        }

        free(list);
        free(nu);
        free(w);
        free(u);
        free(tmp);
        free(stab);

}

int property(int *g, int n, int prop, int *info, int infol) {

        int *perm1= (int*)malloc(n*sizeof(int));
        int *perm2= (int*)malloc(n*sizeof(int));
        int i;
        int ret=0; /* By default we do not accept g */

        if (prop==1) { /* Centralizer. info is the permutation to centralize
                          infol is the length of the perm cycle used in the base. Not used to check property */
                product(g,info,perm1,n);
                product(info,g,perm2,n);
                ret=equal_list(perm1,perm2,n);
        } else if (prop==2) { /* Normalizer */
        } else if (prop==3) { /* Group intersection. */
        } else if (prop==4) { /* SetStabilizer. info is the list of set-points as a characteristic function
                                 infol is the number of points in the original list, not in the characteristic function (which is n) */
                for(i=0; i<n; i++) {
                        if(info[i]) {
                                ret=info[onpoints(i+1,g,n)-1];
                                if(!ret) break;
                        }
                }
        }

        free(perm1);
        free(perm2);
        
        return ret;

}

void generate(int *base, int bl, int *GS, int m, int n, int prop, int *info, int infol,
        int s, int i, int *list, int ll, int **GSK, int *mK, int *mark, int *num) {

        int *g=     (int *)malloc(  n*sizeof(int));
        int *stab=  (int *)malloc(m*n*sizeof(int)), sl;
        int *orbit= (int *)malloc(  n*sizeof(int)), ol;
        int *tmp=   (int *)malloc(  n*sizeof(int));
        int j,gamma;

/*
        FILE *file;
*/

        from_base_image(list, ll, base, bl, GS, m, n, g);

        if (i==bl+1) {
                (*num)++;
/*
        file = fopen("/home/jmm/setwise","a");
        fprintf(file,"generate at level s=%d and i=%d. We have list=",s,i);
        fprint_list(file, list, ll, 0);
        fprintf(file," and permutation ");
        fprint_perm(file, g, n, 1);
        fclose(file);
*/
                if(!isid(g,n) && property(g,n,prop,info,infol)) {
                        /* Append permutation g to SGSK */
                        if (*mK>=m) {
                                *GSK = realloc(*GSK, (*mK+1)*n*sizeof(int));
                        }
                        copy_list(g,(*GSK)+(*mK)*n,n);
                        (*mK)++;
                        *mark=1;
                }
        } else {
                stabilizer(base,i-1,GS,m,n, stab, &sl);
                one_orbit(base[i-1],stab,sl,n, orbit, &ol);
/*
        file = fopen("/home/jmm/setwise","a");
        fprintf(file,"generating over points ");
        fprint_list(file, orbit, ol, 1);
        fclose(file);
*/
                for (j=0; j<ol; j++) {
                        gamma = onpoints(orbit[j], g, n);
                        if(prop==1 && onpoints(list[ll-1],info,n)!=gamma && ll<infol) continue;
/* This needs to be commented out because otherwise this case does not work:
SetStabilizer[{1, 4}, StrongGenSet[{1, 3},
       GenSet[Cycles[{1, 2}], Cycles[{3, 4}], Cycles[{1, 3}, {2, 4}]]] */
                        /* if(prop==4 && !info[gamma-1] && ll<infol) continue; */
                        copy_list(list, tmp, ll);
                        tmp[ll]=gamma;
                        generate(base,bl,GS,m,n, prop,info,infol, s, i+1, tmp, ll+1, GSK, mK, mark, num);
                        if (*mark) break;
                }
        }

        free(g);
        free(stab);
        free(orbit);
        free(tmp);

}

void search(int *base, int bl, int *GS, int m, int n, int prop, int *info, int infol,
        int s, int **GSK, int *mK, int *num) {

        int *stab=  (int *)malloc(m*n*sizeof(int)), sl;
        int *orbit= (int *)malloc(  n*sizeof(int)), ol;
        int *tmp=   (int *)malloc(  n*sizeof(int));
        int *Korbit=(int *)malloc(  n*sizeof(int)), Kol;
        int *list  =(int *)malloc(  n*sizeof(int)), ll;
        int gammas, minKorbit, i, mark;
/*
        FILE *file;

        file = fopen("/home/jmm/setwise","a");
        fprintf(file,"search: s=%d with %d points in base and degree %d\n",s,bl,n);
        fclose(file);
*/
        if (s==bl+1) {
                /* Initialize with the identity group */
                *mK=0;
        } else {
                search(base, bl, GS, m, n, prop,info,infol, s+1, GSK, mK, num);
                stabilizer(base, s-1, GS, m, n, stab, &sl);
                one_orbit(base[s-1],stab,sl,n,orbit,&ol);
/*
        file = fopen("/home/jmm/setwise","a");
        fprintf(file,"branching over points ");
        fprint_list(file, orbit, ol, 1);
        fclose(file);
*/
                /* Note that we do not consider the first point to avoid rechecking permutations. Not in Butler */
                for (i=1; i<ol; i++) {
                        gammas=orbit[i];
                        one_orbit(gammas,*GSK,*mK,n,Korbit,&Kol);
                        sortB(Korbit,tmp, Kol, base, bl);
                        minKorbit=tmp[0];
/*
        file = fopen("/home/jmm/setwise","a");
        fprintf(file,"search: i=%d, gammas=%d, minKorbit=%d\n",i, gammas, minKorbit);
        fclose(file);
*/
        
                        if (minKorbit==gammas) {
                                /* First point in orbit */
                                copy_list(base,list,s-1);
                                list[s-1]=gammas;
                                ll=s;
                                mark=0;
                                generate(base,bl,GS,m,n,prop,info,infol, s,s+1,list,ll, GSK,mK,&mark, num);
                        }
                }
        }

        free(stab);
        free(orbit);
        free(tmp);
        free(Korbit);
        free(list);

}

/* stab_chain. TB, 1 February 2014
 *
 * Computes the stabilizer chain from the SGS given by
 * base (length bl) and GS (m n-permutations)
 * The result is chain a list of pointers to lists 
 * containing the positions of the generators within GS.
 * The lengths of the lists are given by cl, which is a list 
 * with length bl. Observe that we do not save the last elemet 
 * in the chain which should always be empty.
 * We assume that enough space is already allocated for cl 
 * (bl integers), and for chain (bl pointers) 
 * 
 */

void stab_chain(int *base, int bl, int *GS, int m, int n, int **chain, int *cl) {

	int i, j;
	int gsp[m]; 

/* The first GenSet should be the full list.
 */

	cl[0] = m;
	chain[0] = malloc(m*sizeof(int));
	
	j = m;
	while(j--) chain[0][j] = j;


/* For the rest we need to check which perms that fixes the 
   previous base point and is in the previous stabilizer.
 */

	for(i = 1; i<bl; i++) {
		for(j=0; j<cl[i-1]; j++) {
			if (onpoints(base[i-1],&GS[n*(chain[i-1][j])],n) == base[i-1]) {
				gsp[cl[i]]=chain[i-1][j];
				cl[i]++;
			}
		}
		chain[i] = malloc(cl[i]*sizeof(int));
		memmove(chain[i], gsp, cl[i]*sizeof(int));
	} 
}


/* one_orbit_chain. TB, 16 March 2014
 *
 * Does the same thing as one_orbit, but instead of giving the 
 * generating set as a simple list, it is given as a list GS together with
 * a list of positions poslist with length poslistl. 
 *
 */


void one_orbit_chain(int point, int *GS, int *poslist, int poslistl, 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 poslist */
	int newgamma;

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

}


/* one_schreier_orbit_chain. TB, 16 March 2014
 *
 * Does the same thing as one_schreier_orbit, but instead of giving the 
 * generating set as a simple list, it is given as a list GS together with
 * a list of positions poslist with length poslistl. 
 *
 */


void one_schreier_orbit_chain(int point, int *GS, int *poslist, int poslistl,
	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 poslistl */
	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<poslistl; mp++) {
			copy_list(&(GS[poslist[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);

}



/* conjugate_chain. TB, 16 February 2014
 *
 * Conjugates the stabilizer chain with respect to the 
 * permutation p.
 * 
 */

void conjugate_chain(int *base, int bl, int *GS, int m, int n, int *p) {
	int j;
	int tmpperm[n];
	int ip[n];
	
	inverse(p, ip, n);


/* Transform the base */
	for(j = 0; j < bl; j++) {
		base[j] = onpoints(base[j], p, n);
	}

/* Transform the GenSet */
	for(j = 0; j < m; j++) {
		product(&GS[n*j], p, tmpperm, n);
		product(ip, tmpperm, &GS[n*j], n);
	}

}


/* basechange_chain. TB, 12 March 2014
 *
 * Changes the base in the stabilizer chain.
 * 
 */

void basechange_chain(int **base, int *bl, int **GS, int *m, int n, int ***chain, int **cl, int *newbase, int newbl) {

	int orbit[n], ol;
	int more=0;
	int nu[n*n];
	int w[n];
	int g[n];
	int invg[n];
	int g2[n];
	int g3[n];
	int i,j, point;
/*	int GSK[(*m)*n];*/
	int pos;
/*	FILE *file;*/
	
	i=1;

	if(newbl>0){
		one_schreier_orbit((*base)[0], *GS, *m, n, orbit, &ol, nu, w, 1);
		more=position(newbase[0], orbit, ol);
	}
	else{
		more=0;
	}
	if(more>0){
		trace_schreier(newbase[0], nu, w, g, n);
		while(more>0 && i<*bl && i<newbl){
/*			for(j=0; j < (*cl)[i]; j++) {
				memmove(&GSK[n*j], &(*GS)[n*((*chain)[i][j])], n*sizeof(int));
			}*/
			inverse(g, invg, n);
			point=onpoints(newbase[i], invg, n);
/*			one_schreier_orbit((*base)[i], GSK, (*cl)[i], n, orbit, &ol, nu, w, 1);*/
			one_schreier_orbit_chain((*base)[i], *GS, (*chain)[i], (*cl)[i], n, orbit, &ol, nu, w, 1);
			more=position(point, orbit, ol);
			if(more>0){
				trace_schreier(point, nu, w, g2, n);
				product(g2,g,g3,n);
				memmove(g, g3, n*sizeof(int));
			} 
			i=i+1;
		}
	conjugate_chain(*base, *bl, *GS, *m, n, g); 
	}
	for(j=i-1; j < newbl; j++) {
		pos=position(newbase[j], *base, *bl);
		if(pos==0){
			appendbasepoint_chain(base, bl, chain, cl, newbase[j]);
			pos=*bl;
		}
/*	        file = fopen("xpermlog.txt","a");
        	fprintf(file,"i=%d, j=%d, pos=%d \n", i, j, pos);
		fprint_list(file, *base, *bl, 1);
        	fclose(file);*/

		while((pos!=(j+1))&&(pos>1)){
/*		        file = fopen("xpermlog.txt","a");
	        	fprintf(file,"interchange_chain, pos=%d \n", pos);
        		fclose(file);*/

			interchange_chain(base, bl, GS, m, n, chain, cl, --pos);
		}
	}


}

/* appendbasepoint_chain. TB, 12 March 2014
 *
 * Appends a point to the end of the base in the stabilizer chain structure.
 * 
 */

void appendbasepoint_chain(int **base, int *bl, int ***chain, int **cl, int newbasepoint) {

	(*bl)++;
	*base=(int *)realloc(*base, (*bl)*sizeof(int));
	*cl=(int *)realloc(*cl, (*bl)*sizeof(int));
	*chain=(int* *)realloc(*chain, (*bl)*sizeof(int*));
	(*base)[(*bl)-1]=newbasepoint;
	(*cl)[(*bl)-1]=0;
	(*chain)[(*bl)-1]=NULL;

}

/* interchange_chain. TB, 13 March 2014
 *
 * Interchanges the point j and j+1 in the base using the stabilizer chain structure.
 * 
 */

void interchange_chain(int **base, int *bl, int **GS, int *m, int n, int ***chain, int **cl, int j) {
	
	int basej, basejp1;
	int Deltaj[n], Deltajl;
	int nuj[n*n];
	int wj[n];
	int Deltajp1[n], Deltajp1l;
	int orbitjp1jm1l;
	int nujp1[n*n];
	int wjp1[n];
	int i,k;
	int GSK[(*m)*n];
	int LDeltaBarjp1;
	int *T=NULL;
	int Tl=0;
	int *Gamma=NULL;
	int Gammal;
	int Delta[n];
	int Deltal, gamma, p, pos;
	int g1[n];
	int invg1[n];
	int g2[n], g2g1[n];
	int orbitgammaT[n], orbitgammaTl;
	int *newT=NULL;
	int newTl, oldm;
	int *gsp=NULL; 

/*	FILE *file;*/


	basej=(*base)[j-1];
	basejp1=(*base)[j];
	oldm=*m;

/*        file = fopen("xpermlog.txt","a");
       	fprintf(file,"Inside interchange_chain: m=%d, basej=%d, basejp1=%d \n",*m, basej, basejp1);
	fprintf(file,"base=");
	fprint_list(file, *base, *bl, 1);
       	fclose(file);

        file = fopen("xpermlog.txt","a");
	fprintf(file,"cl=");
	fprint_list(file, *cl, *bl, 1);
       	fclose(file);*/


/*	for(i=0; i < (*cl)[j]; i++) {
		memmove(&GSK[n*i], &(*GS)[n*((*chain)[j][i])], n*sizeof(int));
	}
	one_schreier_orbit(basejp1, GSK, (*cl)[j], n, Deltajp1, &Deltajp1l, nujp1, wjp1, 1);
*/

	one_schreier_orbit_chain(basejp1, *GS,(*chain)[j], (*cl)[j], n, Deltajp1, &Deltajp1l, nujp1, wjp1, 1);


/*        file = fopen("xpermlog.txt","a");
	fprintf(file,"Deltajp1=");
	fprint_list(file, Deltajp1, Deltajp1l, 1);
       	fclose(file);*/

	
/*	for(i=0; i < (*cl)[j-1]; i++) {
		memmove(&GSK[n*i], &(*GS)[n*((*chain)[j-1][i])], n*sizeof(int));
	}

	one_orbit(basejp1, GSK, (*cl)[j-1], n, Deltaj, &orbitjp1jm1l);

	one_schreier_orbit(basej, GSK, (*cl)[j-1], n, Deltaj, &Deltajl, nuj, wj, 1);
*/

	one_orbit_chain(basejp1, *GS, (*chain)[j-1], (*cl)[j-1], n, Deltaj, &orbitjp1jm1l);

	one_schreier_orbit_chain(basej, *GS,(*chain)[j-1], (*cl)[j-1], n, Deltaj, &Deltajl, nuj, wj, 1);

/*        file = fopen("xpermlog.txt","a");
	fprintf(file,"Deltaj=");
	fprint_list(file, Deltaj, Deltajl, 1);
       	fclose(file);*/


	LDeltaBarjp1=Deltajp1l*Deltajl/orbitjp1jm1l;	

	if(j+1<*bl){
		T=(int *)malloc(((*cl)[j+1])*n*sizeof(int));
		Tl=(*cl)[j+1];
		for(i=0; i < (*cl)[j+1]; i++) {
			memmove(&T[n*i], &(*GS)[n*((*chain)[j+1][i])], n*sizeof(int));
		}
	}
	
	Gamma=(int *)malloc(Deltajl*sizeof(int));
	Gammal=0;

	sort(Deltaj, Gamma, Deltajl);

	for(i=0; i<Deltajl; i++){
		if((Gamma[i]!=basej)&&(Gamma[i]!=basejp1)){
			Gamma[Gammal++]=Gamma[i];
		}
	}
	Delta[0]=basej;
	Deltal=1;
	
/*        file = fopen("xpermlog.txt","a");
       	fprintf(file,"Delta=");
	fprint_list(file, Delta, Deltal, 1);
       	fclose(file);*/


	while(Deltal < LDeltaBarjp1){

		gamma=Gamma[0];
		trace_schreier(gamma, nuj, wj, g1, n);
		inverse(g1, invg1, n);
		p=onpoints(basejp1, invg1, n);

/*	        file = fopen("xpermlog.txt","a");
       		fprintf(file,"gamma=%d, p=%d\n", gamma, p);
       		fclose(file);*/

		pos=position(p, Deltajp1, Deltajp1l);
		if(pos==0){
/*		        file = fopen("xpermlog.txt","a");
	       		fprintf(file,"p in Deltajp1\n");
       			fclose(file);*/

			one_orbit(gamma, T, Tl, n, orbitgammaT, &orbitgammaTl);
			k=0;
			for(i=0; i<Gammal; i++){
				if(position(Gamma[i],orbitgammaT, orbitgammaTl)==0){
					Gamma[k++]=Gamma[i];
				}
			}
			Gammal=k;
		}
		else{
/*		        file = fopen("xpermlog.txt","a");
	       		fprintf(file,"p not in Deltajp1\n");
       			fclose(file);*/

			trace_schreier(p, nujp1, wjp1, g2, n);
			product(g2,g1,g2g1,n);
			Tl++;
			T=(int *)realloc(T,Tl*n*sizeof(int));
			memmove(&T[n*(Tl-1)], g2g1, n*sizeof(int));
			one_orbit(basej, T, Tl, n, Delta, &Deltal);

			
			k=0;
			for(i=0; i<Gammal; i++){
				if(position(Gamma[i], Delta, Deltal)==0){
					Gamma[k++]=Gamma[i];
				}
			}
			Gammal=k;

		}

		
	}

/*        file = fopen("xpermlog.txt","a");
       	fprintf(file,"After loop: Tl=%d\n", Tl);
	fprint_list(file,T,Tl*n,1);
       	fclose(file);*/

	newT=(int *)malloc(Tl*n*sizeof(int));
	newTl=0;

	for(i=0; i < (*cl)[j-1]; i++) {
		memmove(&GSK[n*i], &(*GS)[n*((*chain)[j-1][i])], n*sizeof(int));
	}

	if(j+1<*bl){
		complement(&T[n*((*cl)[j+1])], Tl-(*cl)[j+1], GSK, (*cl)[j-1], n, newT, &newTl);
	}
	else{
		complement(T, Tl, GSK, (*cl)[j-1], n, newT, &newTl);
	}
	

	(*GS)=(int *)realloc(*GS, (newTl+oldm)*n*sizeof(int));
	memmove(&((*GS)[n*(*m)]), newT, newTl*n*sizeof(int));
	(*m)=newTl+oldm;

	(*base)[j-1]=basejp1;
	(*base)[j]=basej;

/*        file = fopen("xpermlog.txt","a");
       	fprintf(file,"newTl=%d\n", newTl);
	fprint_list(file,newT,newTl*n,1);
       	fclose(file);*/


	for(i=0;i<j;i++){
		(*chain)[i] = realloc((*chain)[i],(newTl+(*cl)[i])*sizeof(int));
		for(k=0;k<newTl;k++){
			(*chain)[i][k+(*cl)[i]]=k+oldm;
		}
		(*cl)[i]=newTl+(*cl)[i];
		
	}

	gsp=malloc(((*cl)[j-1])*sizeof(int));

	(*cl)[j]=0;

	for(i=0; i<(*cl)[j-1]; i++) {
		if (onpoints(basejp1, &((*GS)[n*((*chain)[j-1][i])]),n) == basejp1) {
			gsp[(*cl)[j]]=(*chain)[j-1][i];
			(*cl)[j]++;
		}
	}
	(*chain)[j] = realloc((*chain)[j],((*cl)[j])*sizeof(int));
	memmove((*chain)[j], gsp, (*cl)[j]*sizeof(int));

	free(gsp);

	free(T);
	free(newT);
	free(Gamma);

}