/********************************************************************/
/*                                                                  */
/*  Module        : Cohomology                                      */
/*                                                                  */
/*  Version       : 2.2                                             */
/*  Last revision : 03/30/94 15:02:14                               */
/*                                                                  */
/*  Description :                                                   */
/*     This module is used to compute the first and second coho-    */
/*     mology groups of agiven p-group.                             */
/*                                                                  */
/*  Functions supplied :                                            */
/*     int calc_h1 ( int dimension );                               */
/*                                                                  */
/********************************************************************/

#include "aglobals.h"
#include "fdecla.h"
#include	"pc.h"
# include	"aut.h"
#include	"hgroup.h"
#include	"grpring.h"
#include	"storage.h"

void az1_mat 			_(( int homogeneous ));
PCGRPDESC *p_quotient  	_(( PCGRPDESC *g_desc, int class ));
VEC *ig_op			_(( void ));
VEC *tensor_prod 		_(( VEC rep1[], VEC rep2[], int dim1, int dim2 ));
void calc_h1 			_(( VEC g_op[], int dimension, int *z1_dim, int *h1_dim ));
VEC *triv_mod 			_(( int dim, int numgen ));
void cohomology 		_(( PCGRPDESC *g_desc, int n, VEC g_op[], int dimension, int *zn_dim, int *hn_dim ));
VEC expand_two_cocycle 	_(( VEC z1, PCELEM e1, PCELEM e2 ));

extern PCGRPDESC *group_desc;
extern GRPDSC *h_desc;
extern int start, bperelem;
extern char matrix[YMAX][XMAX];
extern int x_dim, y_dim;
extern VEC absolut, inhom;
extern VEC fsolution[XMAX];
extern int dim, dquad;
extern int adim, adquad;
extern PCELEM *rho;
extern VEC *opmatrix;	
extern HOM *dgroup_auts;

VEC *calc_b1 ( VEC g_op[], int dimension, int *b1_dim )
{
	VEC *b1;
	VEC m_id, temp;
	int i, j, k, dq, xs, ys;

	PUSH_STACK();
	dq = dimension * dimension;
	m_id = CALLOCATE ( dq );
	for ( i = 0; i < dimension; i++ )
		m_id[i*dimension+i] = GPRIME-1;
	for ( k = 0; k < bperelem; k++ ) {
		temp = ALLOCATE ( dq );
		copy_vector ( g_op[k], temp, dq );
		ADD_VECTOR ( m_id, temp, dq );
		for ( i = 0; i < dimension; i++ ) {
			for ( j = 0; j < dimension; j++ )
				matrix[(long)i][(long)j+k*dimension] = temp[i*dimension+j];
		}
	}
	xs = x_dim;
	ys = y_dim;
	x_dim = dimension * bperelem;
	y_dim = dimension;
	absolut = CALLOCATE ( y_dim );
	inhom = ALLOCATE ( x_dim );
	k = GAUSS_ELIMINATE();
	POP_STACK();

	b1 = ARRAY ( k, VEC );
	j = 0;
	for ( i = 0; i < dimension; i++ ) {
		if ( !iszero ( matrix[(long)i], x_dim ) ) {
			b1[j] = ALLOCATE ( x_dim );
			copy_vector ( matrix[(long)i], b1[j++], x_dim );
		}
	}
			
	/* restore old values of x_dim and y_dim */
	x_dim = xs;
	y_dim = ys;
	*b1_dim = k;
	return ( b1 );
}

VEC *calc_z1 ( VEC g_op[], int dimension, int *z1_dim )
{
	int i, j, k, xd, yd;
	int sdim, sdquad;
	VEC *z1;
	
	start = bperelem;

	PUSH_STACK();
	
	/* save old values of dim and dquad */
	sdim = dim;
	sdquad = dquad;
	dim = adim = dimension;
	dquad = adquad = adim * adim;

	/* compute operating matrices for identity */
	opmatrix = ARRAY ( start, VEC );
	for ( j = start; j--; ) {
		opmatrix[j] = g_op[j];
	}

	absolut = CALLOCATE ( adim * NUMREL );
	inhom = ALLOCATE ( adim * bperelem );
	rho = ARRAY ( bperelem, PCELEM );
	
	for ( i = 0; i < bperelem; i++ ) {
		rho[i] = IDENTITY;
		rho[i][i] = 1;
	}

	az1_mat ( TRUE );
	xd = bperelem * adim;
	yd = NUMREL * adim;
	*z1_dim = xd - solve_equations ( xd, yd );
	
	k = 0;
	for ( i = xd; i--; ) {
		if ( fsolution[i] ) {
			copy_vector ( fsolution[i], matrix[(long)k], xd );
			k++;
		}
	}
	POP_STACK();

	z1 = ARRAY ( *z1_dim, VEC );
	for ( i = 0; i < *z1_dim; i++ ) {
		z1[i] = ALLOCATE ( xd );
		copy_vector ( matrix[(long)i], z1[i], xd );
	}
			
	/* restore old values of dim and dquad */
	dim = sdim;
	dquad = sdquad;

	return ( z1 );
}

void calc_h1 ( VEC g_op[], int dimension, int *z1_dim, int *h1_dim )
/* compute largerst i with P_class/P_i elementary abelian and setup
   operating matrices for G/P_class operating on P_class/P_i.
   Setup and solve system of linear equations.
*/
{
	int i, j, k;
	int sdim, sdquad;
	int xs, ys, xd, yd;
	VEC *h1;
	VEC m_id, temp;
	
	start = bperelem;

	PUSH_STACK();
	
	/* save old values of dim and dquad */
	sdim = dim;
	sdquad = dquad;
	dim = adim = dimension;
	dquad = adquad = adim * adim;

	/* compute operating matrices for identity */
	opmatrix = ARRAY ( start, VEC );
	for ( j = start; j--; ) {
		opmatrix[j] = g_op[j];
	}

	absolut = CALLOCATE ( dim * NUMREL );
	inhom = ALLOCATE ( NUMGEN * bperelem );
	rho = ARRAY ( bperelem, PCELEM );
	
	for ( i = 0; i < bperelem; i++ ) {
		rho[i] = IDENTITY;
		rho[i][i] = 1;
	}

/*	for ( j = 0; j < start; j++ ) {
		printf ( "\nmat. no, %d\n", j );
		show_mat ( opmatrix[j] );
	}
*/
	az1_mat ( TRUE );
	xd = bperelem * adim;
	yd = NUMREL * adim;
	solve_equations ( xd, yd );
	
	m_id = CALLOCATE ( dquad );
	for ( i = 0; i < dim; i++ )
		m_id[i*dim+i] = GPRIME-1;
	for ( k = 0; k < bperelem; k++ ) {
		temp = ALLOCATE ( dquad );
		copy_vector ( g_op[k], temp, dquad );
		ADD_VECTOR ( m_id, temp, dquad );
		for ( i = 0; i < dim; i++ ) {
			for ( j = 0; j < dim; j++ )
				matrix[(long)i][(long)j+k*dim] = temp[i*dim+j];
		}
	}

	k = dim;
	for ( i = xd; i--; ) {
		if ( fsolution[i] ) {
			copy_vector ( fsolution[i], matrix[(long)k], xd );
			k++;
		}
	}
	POP_STACK();
	*z1_dim = k - dim;

	xs = x_dim;
	ys = y_dim;
	x_dim = xd;
	y_dim = k;

	*h1_dim = complement ( dim, xd, k );

	h1 = ARRAY ( *h1_dim, VEC );
	for ( i = 0; i < *h1_dim; i++ ) {
		h1[i] = fsolution[i];
	}
			
	/* restore old values of dim and dquad */
	dim = sdim;
	dquad = sdquad;

	x_dim = xs;
	y_dim = ys;
}

	
static int dimension, ig_dim;
static int c2len;
static VEC zindex;
static VEC *t;

VEC expand_two_cocycle ( VEC z1, PCELEM e1, PCELEM e2 )
{
	int i, j, k, len;
	VEC val;
	PCELEM f1, f2;
	
	val = CALLOCATE ( dimension );
	if ( iszero ( e1, bperelem ) || iszero ( e2, bperelem ) )
		return ( val );
	i = 0;
	while ( e1[i] == 0 ) i++;
	j = bperelem - 1;
	while ( e1[j] == 0 ) j--;
	len = j - i + 1;
	if ( len == 1 ) {
		k = IND ( g_invers ( e2 ) );
		copy_vector ( z1+i*dimension*ig_dim+(k-1)*dimension, val, dimension );
	}
	else {
		PUSH_STACK();
		f1 = IDENTITY;
		copy_vector ( e1, f1, i+1 );
		f2 = IDENTITY;
		copy_vector ( e1+i+1, f2+i+1, bperelem - i - 1 );
		copy_vector ( expand_two_cocycle ( z1, f2, e2 ), val, dimension );
		ADD_VECTOR ( expand_two_cocycle ( z1, f1, monom_mul ( f2, e2 ) ), val, dimension );
		SUBB_VECTOR ( expand_two_cocycle ( z1, f1, f2 ), val, dimension );
		POP_STACK();
	}
	return ( val );
}
		
VEC corr_one_cocycle ( VEC fs, int dim )
{
	int d = dim * (GCARD - 1);
	int i, k;
	PCELEM el2;
	VEC z1;
	
	z1 = ALLOCATE ( bperelem * d );
	
	PUSH_STACK();
	el2 = IDENTITY;
	while ( inc_count ( el2, bperelem ) ) {
		for ( i = 0; i < bperelem; i++ ) {
			k = IND ( g_invers ( el2 ) );
			copy_vector ( fs+((k-1)*bperelem+i)*dim, z1+i*d+(k-1)*dim, dim );
		}
	}
	POP_STACK();
	return ( z1 );
}	
	
VEC corr_two_cocycle ( VEC z1, VEC alpha )
{
	int i, j;
	PCELEM el1, el2;
	VEC value;
	VEC fs;
	
	fs = ALLOCATE ( dimension * c2len );
	
	PUSH_STACK();
	el1 = IDENTITY;
	el2 = IDENTITY;
	j = 0;
	while ( inc_count ( el2, bperelem ) ) {
		PUSH_STACK();
		for ( i = 0; i < bperelem; i++ ) {
			zero_vector ( el1, bperelem );
			el1[i] = 1;
			if ( alpha == NULL )
				value = expand_two_cocycle ( z1, el1, el2  );
			else	
				value = expand_two_cocycle ( z1, image ( alpha, el1 ), image ( alpha, el2 ) );
			copy_vector ( value, fs + j*dimension, dimension );
			j++;
		}
		POP_STACK();
	}
	POP_STACK();
	return ( fs );
}

VEC factor_set ( PCGRPDESC *g_desc )
{
	PCGRPDESC *old_pc_group;
	PCELEM el1, el2, r1, r2;
	int i, j, start, dim, entries;
	VEC fs;

	old_pc_group = group_desc;
	set_main_group ( g_desc );
	
	start = EXP_P_LCS[EXP_P_CLASS].i_start;
	dim = GNUMGEN - start;
	entries = start;
	for ( i = 0; i < start; i++ )
		entries *= GPRIME;
	entries -= start;
	fs = ALLOCATE ( dim * entries );

	PUSH_STACK();
	el1 = IDENTITY;
	el2 = IDENTITY;
	j = 0;
	while ( inc_count ( el2, start ) ) {
		PUSH_STACK();
		for ( i = 0; i < start; i++ ) {
			zero_vector ( el1, bperelem );
			el1[i] = 1;
			r1 = monom_mul ( el1, el2 );
			r2 = IDENTITY;
			copy_vector ( r1, r2, start );
			r1 = monom_mul ( g_invers ( r2 ) , r1 );
			copy_vector ( r1+start, fs + j*dim, dim );
			j++;
		}
		POP_STACK();
	}
	POP_STACK();
	set_main_group ( old_pc_group );
	return ( fs );
}

VEC *calc_c2 ( VEC z1[], int z1_dim )
{
	int i, j, k, entries;
	PCELEM el1, el2;
	VEC value;
	VEC *fs;
	
	entries = bperelem;
	for ( i = 0; i < bperelem; i++ )
		entries *= GPRIME;
	entries -= bperelem;
	
	c2len = entries;
	fs = ARRAY ( z1_dim, VEC );
	for ( i = 0; i < z1_dim; i++ )
		fs[i] = ALLOCATE ( dimension * entries );
	
	for ( k = 0; k < z1_dim; k++ ) {
		PUSH_STACK();
		el1 = IDENTITY;
		el2 = IDENTITY;
		j = 0;
		while ( inc_count ( el2, bperelem ) ) {
			PUSH_STACK();
			for ( i = 0; i < bperelem; i++ ) {
				zero_vector ( el1, bperelem );
				el1[i] = 1;
				value = expand_two_cocycle ( z1[k], el1, el2 );
				copy_vector ( value, fs[k] + j*dimension, dimension );
				j++;
			}
			POP_STACK();
		}
		POP_STACK();
	}
	return ( fs );
}

void get_trafo_mat ( VEC *b2l, int bdim )
{
	int i, j;
	int xd, yd;
	
	xd = bdim;
	yd = c2len;

	PUSH_STACK();
	for ( i = 0; i < yd; i++ )
		for ( j = 0; j < xd; j++ )
			matrix[(long)i][(long)j] = b2l[j][i];
	gauss_p_eliminate ( xd, yd );
	POP_STACK();
	
	t = ARRAY ( yd, VEC );
	zindex = CALLOCATE ( yd );
	for ( i = 0; i < yd; i++ ) {
		t[i] = ALLOCATE ( yd );
		copy_vector ( matrix[(long)i]+(long)xd, t[i], yd );
		zindex[i] = ( iszero ( matrix[(long)i], xd ) );
	}
}

void	addlist ( VEC *orbit, int *offset, VEC fs, int bdim )
{
	int i, j, k, isnew;
	int xd, yd;
	register int val;
	VEC tt;
	
	PUSH_STACK();
	xd = bdim;
	yd = c2len;
	
	absolut = ALLOCATE ( yd );
	inhom = ALLOCATE ( xd );

	isnew = TRUE;

/*	for ( k = 0; k < *offset; k++ ) {
		copy_vector ( fs, absolut, c2len );
		SUBB_VECTOR ( orbit[k], absolut, c2len );
		for ( i = 0; i < yd; i++ )
			for ( j = 0; j < xd; j++ )
				matrix[(long)i][(long)j] = b2l[j][i];
		
		isnew = ( solve_equations ( xd, yd ) == -1 );
		if ( !isnew )
			break;
	} */
	
	for ( k = 0; k < *offset; k++ ) {
		isnew = FALSE;
		copy_vector ( fs, absolut, c2len );
		SUBB_VECTOR ( orbit[k], absolut, c2len );
		for ( i = 0; i < yd; i++ ) {
			if ( zindex[i] ) {
				tt = t[i];
				val = 0;
				for ( j = 0; j < yd; j++ )
					val += (tt[j]*absolut[j]);
/*					val = ADD ( val, MUL ( t[i][j], absolut[j] ) ); */
				if ( (isnew = ( (val % GPRIME)!= 0 )) ) break;
			}
			
		}
		if ( !isnew ) break; 
	} 
	
	POP_STACK();
	
	if ( isnew ) {
		orbit[*offset] = ALLOCATE ( c2len );
		copy_vector ( fs, orbit[*offset], c2len );
		++(*offset);
	}
}

void two_coboundaries ( PCGRPDESC *g_desc, VEC g_op[], int dim )
{
	VEC *aug;
	VEC *ng_op;
	VEC *b1, *b2, *orbit, z1fs, alpha, fs;
	int b1_dim, orbit_len, i, auts;
	PCGRPDESC *old_pc_group, *gg;
	GRPDSC *old_p_group;

	fs = factor_set ( g_desc );

	old_pc_group = group_desc;
	old_p_group = h_desc;
	gg = p_quotient ( g_desc, g_desc->exp_p_class - 1);
	set_main_group ( gg );
	set_h_group (  conv_rel ( gg ) );

	/* compute B^1(g, Hom(IG,A)) */
	aug = ig_op();
	ig_dim = GCARD - 1;
	ng_op = tensor_prod ( aug, g_op, ig_dim, dim );
	dimension = dim * ig_dim;
	b1 = calc_b1 ( ng_op, dimension, &b1_dim );

	/* compute elements of B^2(G,A) */
	dimension = 1;
	b2 = calc_c2 ( b1, b1_dim );

	get_trafo_mat ( b2, b1_dim );
	z1fs = corr_one_cocycle ( fs, dimension );
	dgroup_auts = automorphisms ( gg, 0 );
	dgroup_auts = generate_automorphism_group ( dgroup_auts, FALSE );

	auts = dgroup_auts->aut_gens_dim[1];
	orbit = ARRAY ( auts, VEC );
	orbit_len = 0;
	
	for ( i = 0; i < auts; i++ ) {
		alpha = dgroup_auts->aut_gens[1][i];
		fs = corr_two_cocycle ( z1fs, alpha );
		addlist ( orbit, &orbit_len, fs, b1_dim );
	}
	
	printf ( "length of orbit: %d\n", orbit_len );
/*	for ( i = 0; i < orbit_len; i++ )
		write_vector ( orbit[i], c2len ); */
	set_main_group ( old_pc_group );
	set_h_group ( old_p_group );
}

void cohomology ( PCGRPDESC *g_desc, int n, VEC g_op[], int dimension, int *zn_dim, int *hn_dim )
{
	VEC *aug;
	VEC *ng_op;
	int d_aug, dim;
	int i;
	PCGRPDESC *old_pc_group;
	GRPDSC *old_p_group;

	old_pc_group = group_desc;
	old_p_group = h_desc;
	set_main_group ( g_desc );
	set_h_group (  conv_rel ( g_desc ) );

	if ( n == 1 )
		calc_h1 ( g_op, dimension, zn_dim, hn_dim );
	else {
		aug = ig_op();
		d_aug = GCARD - 1;
		dim = dimension;
		ng_op = g_op;
		for ( i = 2; i <= n; i++ ) {
			ng_op = tensor_prod ( aug, ng_op, d_aug, dim );
			dim *= d_aug;
			printf ( "dimension of g-module : %d\n", dim );
		}
/*		for ( j = 0; j < num_gen; j++ ) {
			showmat ( ng_op[j], dim, dim );
			puts ( "\n" );
		} */
		calc_h1 ( ng_op, dim, zn_dim, hn_dim );
	}
	set_main_group ( old_pc_group );
	set_h_group ( old_p_group );
}

void calc_cohomology ( int n, PCGRPDESC *g )
{
	int dim2;
	int zd, hd;
	VEC *g_op;
	PCGRPDESC *gg;

	PUSH_STACK();
	dim2 = g->exp_p_lcs[g->exp_p_class].i_dim;
	gg = p_quotient ( g, g->exp_p_class - 1);
	g_op = triv_mod ( dim2, gg->num_gen );
	
	cohomology ( gg, n, g_op, dim2, &zd, &hd ); 
	printf ( "dimension of Z%1d : %d\n", n, zd );
	printf ( "dimension of H%1d : %d\n", n, hd );
	POP_STACK();
}

void calc_extorbit ( PCGRPDESC *g )
{
	int dim2;
	VEC *g_op;

	PUSH_STACK();
	dim2 = g->exp_p_lcs[g->exp_p_class].i_dim;
	g_op = triv_mod ( dim2, g->num_gen );
	
	two_coboundaries ( g, g_op, dim2 );
	POP_STACK();
}

/* end of modulo cohomology */