/* file: matrix.c */ /* AK 091086 */
#include "def.h"
#include "macro.h"

#ifdef MATRIXTRUE
static struct matrix * callocmatrix();
static INT scan_matrix_co();

static transform_matrix(a,f,b) OP a,b; INT (*f)();
{
	INT e = 0L;
	INT i,j,erg = OK;
	if (a==b)
		{	
		OP c = callocobject();
		*c = *a;
		C_O_K(b,EMPTY);
		e = 1L;
		a = c;
		}
	m_ilih_m(S_M_LI(a),S_M_HI(a),b);
	for (i=0L;i<S_M_HI(b);i++)
		for (j=0L;j<S_M_LI(b);j++)
			erg += (*f)(S_M_IJ(a,i,j),S_M_IJ(b,i,j));
	if (e==1L)
		erg += SYMfreeall(a);
	ENDR("internal function:transform_matrix");
}

INT cast_apply_matrix(a) OP a;
/* AK 270295 */
{
	INT i,ml,erg = OK,j;
	OP b;
	if (S_O_K(a) == MATRIX)
		return OK;
	else if (S_O_K(a) == KOSTKA)
		return OK;
	else if (S_O_K(a) == VECTOR)
		{
		ml = 0;
		for (i=(INT)0;i<S_V_LI(a);i++)
			{
			if (not VECTORP(S_V_I(a,i)))
				goto came;
			if (S_V_LI(S_V_I(a,i)) > ml)
				ml = S_V_LI(S_V_I(a,i));
			}
			
		/* now vector of vector == we can cast */
		b = callocobject();
		*b = *a;
		C_O_K(a,EMPTY);
		erg += m_ilih_m(ml,S_V_LI(b),a);
		for (i=0;i<S_M_HI(a);i++)
		for (j=0;j<S_V_LI(S_V_I(b,i));j++)
			{
			* S_M_IJ(a,i,j) = * S_V_I(S_V_I(b,i),j);
			C_O_K(S_V_I(S_V_I(b,i),j),EMPTY);
			}
		erg += SYMfreeall(b);
		goto camb;
		}
came:
	printobjectkind(a);
	erg += error("cast_apply_matrix: transfer not possible");
camb:
	ENDR("cast_apply_matrix");
}

INT mem_size_matrix(a) OP a;
/* AK 150295 */
{
	INT erg = 0;
        if (a == NULL)
                return 0;
        if (not MATRIXP(a))
                WTO("mem_size_matrix",a);
        erg += sizeof(struct object);
        erg += sizeof(struct matrix);
        erg += mem_size(S_M_H(a));
        erg += mem_size(S_M_L(a));
        erg += S_M_HI(a) * S_M_LI(a) * sizeof(struct object);
        return erg;
}

#endif /* MATRIXTRUE */
#ifdef MATRIXTRUE
INT mod_matrix(a,b,c) OP a,b,c;
/* AK 300793 */
{
	INT erg = OK;
	INT i,j;
	if (S_O_K(b) != INTEGER)
		return ERROR;
	if (not MATRIXP(a))
		return ERROR;
	if (a == c)
		return ERROR;
	if (b == c)
		return ERROR;

	erg += m_ilih_m(S_M_LI(a),S_M_HI(a),c);
	for (i=0L;i<S_M_HI(a);i++)
	for (j=0L;j<S_M_LI(a);j++)
		erg += mod(S_M_IJ(a,i,j),b,S_M_IJ(c,i,j));
	ENDR("mod_matrix");
}

INT matrixp(a) OP a;
/* AK 280193 */
{
	if (S_O_K(a) == MATRIX) 	return TRUE;
	if (S_O_K(a) == KRANZTYPUS) 	return TRUE;
	if (S_O_K(a) == INTEGERMATRIX) 	return TRUE;
	if (S_O_K(a) == KRANZTYPUS) 	return TRUE;
	return FALSE;
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT addinvers_matrix(a,b) OP a,b;
/* AK 240293 */
{
	return transform_matrix(a,addinvers,b);
}

INT absolute_matrix(a,b) OP a,b;
/* AK 240293 */
{
	return transform_matrix(a,absolute,b);
}

INT sum_matrix(a,b) OP a,b;
{
	INT i;
	INT erg = OK;
	OP z = S_M_S(a);
	erg += copy(S_M_IJ(a,S_M_HI(a)-1,S_M_LI(a)-1),b);
	for (i=S_M_HI(a)*S_M_LI(a);i>1;i--)
		erg += add_apply(z++, b);
	ENDR("sum_matrix");
}
	

INT nullp_matrix(a) OP a;
/* AK 110992 */
{
	INT i,j;
	for (i=0L;i<S_M_HI(a);i++)
		for (j=0L;j<S_M_LI(a);j++)
			if (not nullp(S_M_IJ(a,i,j))) return FALSE;
	return TRUE;
}


INT einsp_matrix(a) OP a;
/* AK 110992 */
{
	INT i,j;
	if (S_M_HI(a) != S_M_LI(a)) return FALSE;
	for (i=0L;i<S_M_HI(a);i++)
		for (j=0L;j<S_M_HI(a);j++)
			if (i==j) 
				{
				if (not einsp(S_M_IJ(a,i,j))) return FALSE;
				}
			else	{
				if (not nullp(S_M_IJ(a,i,j))) return FALSE;
				}
	return TRUE;
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT delete_row_matrix(a,index,b) INT index; OP a,b;
/* AK 270789 */ /* AK 111289 V1.1 */ /* AK 110691 V1.2 */ 
/* AK 070891 V1.3 */
{
	INT i,j;
	INT erg = OK;

	if (index >= S_M_HI(a) ) {
		erg += 	error("delete_row_matrix: index to big");
		goto endr_ende;
	}
	if (index < 0L){
		erg += error("delete_row_matrix: index < 0");
		goto endr_ende;
	}

	if (a==b) {
		OP c = callocobject();
		*c = *b; 
		C_O_K(b,EMPTY);
		erg += delete_row_matrix(c,index,b);
		erg += SYMfreeall(c);
		goto endr_ende;
	}

	erg += m_ilih_m(S_M_LI(a),S_M_HI(a)-1L,b);
	for (i=0;i<index;i++)
		for (j=0;j<S_M_LI(a);j++)
			erg += copy(S_M_IJ(a,i,j),S_M_IJ(b,i,j));
	for (i=index+1;i<S_M_HI(a);i++)
		for (j=0;j<S_M_LI(a);j++)
			erg += copy(S_M_IJ(a,i,j),S_M_IJ(b,i-1L,j));
	ENDR("delete_row_matrix");
}
#endif /* MATRIXTRUE */


#ifdef MATRIXTRUE
INT delete_column_matrix(a,index,b) INT index; OP a,b;
/* AK 270789 */ /* AK 310790 V1.1 */ /* AK 110691 V1.2 */ 
/* AK 070891 V1.3 */
{
	INT i,j;
	INT erg = OK;

	if (index >= S_M_LI(a) ) {
		erg += error("delete_column_matrix: index to big");
		goto endr_ende;
	}
	if (index < 0) {
		erg += error("delete_column_matrix: index < 0");
		goto endr_ende;
	}

	if (a==b) {
		OP c = callocobject();
		*c = *b; 
		C_O_K(b,EMPTY);
		erg += delete_column_matrix(c,index,b);
		erg += SYMfreeall(c);
		goto endr_ende;
	}

	erg += m_ilih_m(S_M_LI(a)-1L,S_M_HI(a),b);
	for (j=0;j<index;j++)
		for (i=0;i<S_M_HI(a);i++)
			erg += copy(S_M_IJ(a,i,j),S_M_IJ(b,i,j));
	for (j=index+1;j<S_M_LI(a);j++)
		for (i=0;i<S_M_HI(a);i++)
			erg += copy(S_M_IJ(a,i,j),S_M_IJ(b,i,j-1));
	ENDR("delete_column_matrix");
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT det_matrix(a,b) OP a,b; 
/* AK 151289 V1.1 */ /* AK 110691 V1.2 */
/* hier wird zum konkreten algorithmus geschaltet */
/* AK 210891 V1.3 */
{ 
	if (not quadraticp(a))
		{ 
		error("det_matrix: not quadratic matrix"); 
		return(ERROR); 
		}
	return det_mat_tri(a,b); 
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT det_mat_tri(a,res) OP a,res;
/* AK Fri Jan 20 12:29:58 MEZ 1989
algorithmus zur berechnung der determinante einer matrix a
mittels triangulation
ergebnis in res
algorithmus 41 in CACM */
/* verbessert 1963 */
/* immer noch fehler bei zeilenvertauschung */
/* AK 310790 V1.1 */ /* AK 110691 V1.2 */
/* AK 210891 V1.3 */
{
	INT r,i,j,y,count,sign=1L,n;
	INT erg = OK;
	OP  b,  factor,temp;

	n = S_M_LI(a);
	b = callocobject(); 

	erg += m_i_i(1L,res);
	factor = callocobject(); 
	temp = callocobject();

	erg += copy(a,b);

	for (r=1;r < n;r++)
	{
		count = r-1;
		if (not nullp(S_M_IJ(b,r-1L,r-1))) goto det_mat_tri_resume;
det_mat_tri_zerocheck:
		if (count < (n-1)) count++;
		else goto det_mat_tri_zero;

		if (not nullp(S_M_IJ(b,count,r-1))) {
		for (y=r; y<=n; y++)
			{
			erg += swap(S_M_IJ(b,count,y-1),S_M_IJ(b,r-1L,y-1));
			}
		sign = -sign; goto det_mat_tri_resume;
		}


		goto det_mat_tri_zerocheck;
det_mat_tri_zero:
		erg += m_i_i((INT)0,res); 
		goto det_mat_tri_return;
det_mat_tri_resume:
		for (i=r+1; i<=n; i++)
		{
				
			erg += div3(S_M_IJ(b,i-1L,r-1),S_M_IJ(b,r-1L,r-1),factor);
			for (j=r+1;j<=n;j++)
			{
				erg += mult(factor,S_M_IJ(b,r-1L,j-1),temp);
				erg += addinvers_apply(temp);
				erg += add_apply(temp,S_M_IJ(b,i-1L,j-1L));
			}
		}
	}

	for (i=1;i<=n;i++)
		erg += mult_apply(S_M_IJ(b,i-1L,i-1),res);


	if (sign == -1L) 
		erg += addinvers_apply(res);
	
det_mat_tri_return:
	SYMfreeall(temp); 
	SYMfreeall(factor); 
	SYMfreeall(b);
	ENDR("det_mat_tri");
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT m_ilih_nm(l,h,m) OP m; INT l,h;
/* mit 0 vorbesetzen */
/* make_intlength_intheight_null_matrix */
{
	INT i,erg = OK;
	OP z;
	if (l < 0) { erg += error("m_ilih_nm:l<0");goto endr_ende; }
	if (h < 0) { erg += error("m_ilih_nm:h<0");goto endr_ende; }
	erg += m_ilih_m(l,h,m);
	for (z=S_M_S(m), i=S_M_HI(m) * S_M_LI(m); i>0L; i--,z++)
		M_I_I(0L,z);
	ENDR("m_ilih_nm");
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT m_lh_nm(l,h,m) OP l,h,m;
/* AK 110691 V1.2 */ /* mit 0 vorbesetzen */
/* make_length_height_null_matrix */
/* AK 210891 V1.3 */
{
	INT i,erg = OK;
	OP z;
	CTO(INTEGER,"m_lh_nm",l);
	CTO(INTEGER,"m_lh_nm",h);
	erg += m_lh_m(l,h,m);
	for (z=S_M_S(m), i=S_M_HI(m) * S_M_LI(m); i>0L; i--,z++)
		M_I_I(0L,z);
	ENDR("m_lh_nm");
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT b_lh_nm(l,h,m) OP l,h,m;
/* AK 110691 V1.2 */
/* mit 0 vorbesetzen */
/* build_length_height_null_matrix */
/* AK 210891 V1.3 */
{
	INT i,erg = OK;
	OP z;
	CTO(INTEGER,"b_lh_nm",l);
	CTO(INTEGER,"b_lh_nm",h);
	erg += b_lh_m(l,h,m);
	for (z=S_M_S(m), i=S_M_HI(m) * S_M_LI(m); i>0L; i--,z++)
		M_I_I(0L,z);
	ENDR("b_lh_nm");
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT b_lh_m(l,h,m) OP l,h,m;
/* build_length_height_matrix */
/* height und length werden nicht kopiert */
/* AK 250590 V1.1 */ /* AK 110691 V1.2 */
/* AK 210891 V1.3 */
{
	OP s;
	INT i;
	s = (OP) SYM_malloc(S_I_I(l)*S_I_I(h)*sizeof(struct object));
	for (i=0L;i<S_I_I(l)*S_I_I(h); i++) C_O_K(s+i,EMPTY);
	return b_lhs_m(l,h,s,m);
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT m_lh_m(len,height,matrix) OP height, len, matrix;
/* make_length_height_matrix */
/* height und length werden kopiert */
/* AK 041286 */ /* AK 070789 V1.0 */ /* AK 310790 V1.1 */
/* AK 260291 V1.2 */ /* AK 210891 V1.3 */
{
	INT erg = OK;
	CTO(INTEGER,"m_lh_m",len);
	CTO(INTEGER,"m_lh_m",height);
	if (S_I_I(len) < 0) { erg += error("m_lh_m:negative len");goto endr_ende; }
	if (S_I_I(height) < 0) { erg += error("m_lh_m:negative height");goto endr_ende; }

	erg += b_lhs_m(callocobject(), callocobject(),
	    (struct object *) SYM_calloc((int)
	    ((S_I_I(height))*(S_I_I(len))),
	    sizeof(struct object)),
	    matrix);
	COPY_INTEGER(len,S_M_L(matrix)); 
	COPY_INTEGER(height,S_M_H(matrix));
	ENDR("m_lh_m");
}


INT b_lhs_m(len,height,self,res) OP len, height, self, res;
/* AK 060789 V1.0 */ /* AK 081289 V1.1 */ /* AK 110691 V1.2 */
/* AK 210891 V1.3 */
{
	OBJECTSELF d;

	d.ob_matrix = callocmatrix();
	b_ks_o(MATRIX, d, res);
	C_M_L(res,len); 
	C_M_H(res,height);
	C_M_S(res,self);
	return(OK);
}
#endif /* MATRIXTRUE */


#ifdef MATRIXTRUE
INT m_ilih_m(len,height,matrix) INT height, len; OP matrix;
/* AK 090988 */ /* AK 070789 V1.0 */ /* AK 310790 V1.1 */
/* AK 210891 V1.3 */
{
	INT erg,err;
ma:	erg = OK;
	if (len < 0) {erg += error("m_ilih_m:len < 0"); goto endr_ende;}
	if (height < 0) {erg += error("m_ilih_m:height < 0"); goto endr_ende;}
	erg += b_lhs_m(	callocobject(),
	    callocobject(),
	    (struct object *) SYM_calloc( (int) (height*len),
	    sizeof(struct object)),
	    matrix);

	if (S_M_S(matrix) == NULL) 
		{
		err = error("m_ilih_m:self == NULL ");
		if (err==ERROR_EXPLAIN) {
                        fprintf(stderr,"I wanted a %ld  x %ld matrix", 
				len,height); }
		if (err==ERROR_RETRY) 
			goto ma;
		}
	M_I_I(len,S_M_L(matrix));
	M_I_I(height,S_M_H(matrix));
	ENDR("m_ilih_m");
}


INT quadraticp(mat) OP mat;
/* AK 060789 V1.0 */ /* AK 310790 V1.1 */
/* AK 210891 V1.3 */
{
	return(S_M_LI(mat) == S_M_HI(mat));
}
#endif /* MATRIXTRUE */


#ifdef MATRIXTRUE
#ifdef PERMTRUE
static INT det050995(mat,perm,c) OP mat,perm,c;
/* AK270588 */
/* brechnet aus Matrix (a_ij) und der Permutation p_1,..,p_n
den wert p_1,p_2 * p_3,p_4 * .. */
{
	INT i,erg = OK;

	if (neq(S_M_L(mat),S_P_L(perm)))
		return error("det050995:wrong lengthes");

	erg += copy(S_M_IJ(mat,S_P_II(perm,0L)-1,S_P_II(perm,1L)-1L),c);

	for (i=2L;i<S_P_LI(perm);i+=2)
		{
		erg += mult_apply(S_M_IJ(mat,S_P_II(perm,i)-1L,S_P_II(perm,i+1)-1L),c);
		}

	ENDR("internal routine:det050995");
}
INT det270588(mat,perm,c) OP mat,perm,c;
/* AK270588 */
/* brechnet aus Matrix (a_ij) und der Permutation p_1,..,p_n
den wert a_1,p_1 * a_2,p_2 * .. a_n,p_n */
/* damit kannn z.b. die determinante berechnet werden */
/* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 180691 V1.2 */
/* AK 210891 V1.3 */
{
	INT i,erg = OK;
	if (not quadraticp(mat)) 
		return error("det270588:not quadratic");

	if (neq(S_M_L(mat),S_P_L(perm)))
		return error("det270588:wrong lengthes");

	erg += copy(S_M_IJ(mat,0L,S_P_II(perm,0L)-1L),c);

	for (i=1L;i<S_P_LI(perm);i++)
		{
		erg += mult_apply(S_M_IJ(mat,i,S_P_II(perm,i)-1L),c);
		}

	return erg;
}
#endif /*PERMTRUE */
#endif /* MATRIXTRUE */


#ifdef MATRIXTRUE
#ifdef PERMTRUE
#ifdef CHARTRUE
INT det_mat_imm(mat,erg) OP mat,erg; 
	{ 
	return det_imm_matrix(mat,erg); 
	}

INT det_imm_matrix(mat,b) OP mat,b;
/* AK 270588 neue version mit immanente */
/* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 180691 V1.2 */
/* AK 210891 V1.3 */
{
	OP part;
	INT erg = OK; /* AK 301292 */

	if (not quadraticp(mat)) 
		return error("det_imm_matrix:not quadratic natrir");

	if (S_M_HI(mat) == 1L)  /* AK 301292 */
		{
		erg += copy(S_M_IJ(mat,0L,0L),b);
		goto endr_ende;
		}

	part = callocobject();
	erg += last_partition(S_M_H(mat),part); /* = 1,1,1,1,1,..,1 */
	erg += immanente_matrix(mat,part,b);
	/* der zugehoerige Charakter ist das signum */
	erg += SYMfreeall(part);
	ENDR("det_imm_matrix");
}
#endif /* CHARTRUE */
#endif /* PERMTRUE */
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
#ifdef PERMTRUE
#ifdef CHARTRUE
static INT co_050995(a) OP a;
{
	INT i;
	for (i=0; i< S_P_LI(a) ; i+= 2)
		{
		if (S_P_II(a,i) > S_P_II(a,i+1)) return FALSE;
		if (i > 0)
			if (S_P_II(a,i) < S_P_II(a,i-2)) return FALSE;
		}
	return TRUE;
}

INT pfaffian_matrix(mat,res) OP mat,res;
/* berechnet pfaffian */
/* AK 050995 */
{
	OP perm,zwerg,zzerg;
	INT erg = OK;

	if (not quadraticp(mat)) 
		return error("pffafian:not quadratic matrix");
	if (not evenp(S_M_H(mat))) 
		return error("pffafian:size of matrix not even");

	perm = callocobject(); 
	zwerg = callocobject(); 
	zzerg = callocobject();
	erg += first_permutation(S_M_H(mat),perm);

	erg += det050995(mat,perm,res);
	erg += signum(perm,zwerg);
	erg += mult_apply(zwerg,res);

	while (next(perm,perm))
	{
		if (co_050995(perm) == TRUE)
		{
		erg += det050995(mat,perm,zwerg);
		erg += signum(perm,zzerg);
		erg += mult_apply(zwerg,zzerg);
		erg += add_apply(zzerg,res);
		}
	};

	erg += SYMfreeall(perm);
	erg += SYMfreeall(zwerg);
	erg += SYMfreeall(zzerg);
	ENDR("pfaffian_matrix");
}
INT immanente_matrix(mat,part,erg) OP mat,part,erg;
/* berechnet immanente */
/* AK270588 */ /* AK 060789 V1.0 */ /* AK 090790 V1.1 */ /* AK 180691 V1.2 */
/* AK 210891 V1.3 */
{
	OP perm,nextperm,zwerg,zzerg;

	if (not quadraticp(mat)) 
		return error("immanente:not quadratic matrix");

	perm = callocobject(); 
	zwerg = callocobject(); 
	zzerg = callocobject();
	nextperm = callocobject();
	first_permutation(S_M_H(mat),perm);

	det270588(mat,perm,erg);
	charvalue(part,perm,zwerg,NULL);
	mult_apply(zwerg,erg);

	while (next(perm,nextperm))
	{
		det270588(mat,nextperm,zwerg);
		charvalue(part,nextperm,zzerg,NULL);
		mult_apply(zwerg,zzerg);
		add_apply(zzerg,erg);
		copy(nextperm,perm);
	};

	SYMfreeall(perm);
	SYMfreeall(nextperm);
	SYMfreeall(zwerg);
	SYMfreeall(zzerg);
	return(OK);
}
#endif  /* CHARTRUE */
#endif  /* PERMTRUE */
#endif  /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT inc_matrix(a) OP a;
/* 250488 */ /* AK 060789 V1.0 *//* AK 130790 V1.1 */ /* AK 180691 V1.2 */
/* AK 210891 V1.3 */
{
	OP l=callocobject(),h=callocobject();
	OP b=callocobject(); /* die neue matrix */
	INT i,j;

	COPY_INTEGER(S_M_H(a),h);
	INC_INTEGER(h); 
	COPY_INTEGER(S_M_L(a),l);
	INC_INTEGER(l);
	b_lh_m(l,h,b);C_O_K(b,S_O_K(a));

	for (i=0L;i<S_M_HI(a);i++)
		for (j=0L;j<S_M_LI(a);j++)
			*(S_M_IJ(b,i,j)) =  *(S_M_IJ(a,i,j));
	for (i=0L;i<S_M_HI(b);i++) C_O_K(S_M_IJ(b,i,S_M_LI(a)),EMPTY);
	for (j=0L;j<S_M_LI(b);j++) C_O_K(S_M_IJ(b,S_M_HI(a),j),EMPTY);

	SYM_free(S_M_S(a)); 
	SYMfreeall(S_M_H(a)); 
	SYMfreeall(S_M_L(a));
	SYM_free(S_O_S(a).ob_matrix);

	*a = *b;
	C_O_K(b,EMPTY);
	SYMfreeall(b); 
	return OK;
}
#endif  /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT invers_matrix(a,b) OP a,b;
/* AK 290388 nach stoer (dietmar) */
/* umgewandelt aus pascal */
/* AK 060789 V1.0 */ /* AK 181289 V1.1 */ /* AK 150591 V1.2 */
/* AK 220791 V1.3 */
{
	INT erg=OK;
	INT i,j,k,r;
	/* r ist die selectierte spalte */
	/* r = 0 ... n */
	INT n=S_M_LI(a)-1L;
	INT singulaer = FALSE;
	OP p=callocobject();
	OP hr = callocobject();
	OP hs = callocobject();
	OP hv = callocobject();

	if (not quadraticp(a))
		 return error("invers_matrix:not quadratic"); 

	erg += m_il_v(n+1L,p);

	for(j=0L;j<=n;j++) M_I_I(j,S_V_I(p,j));

	j= -1L;
	erg += copy(a,b);
	while ((j++ <n) && (! singulaer))
	{
		/*pivotsuche*/
		for(r=j;r<=n;r++)
			if (not nullp(S_M_IJ(b,r,j))) goto im290388;

im290388:
		if (r == n+1L)/* nur nullen in der spalte j */ singulaer = TRUE;
		else {
			/*zeilentausch*/
			if (r>j){
				for (k=0L;k<=n;k++) 
					erg+=swap(S_M_IJ(b,j,k),S_M_IJ(b,r,k));
				erg+=swap(S_V_I(p,j),S_V_I(p,r));
			};
			/*transformation*/
			erg+=invers(S_M_IJ(b,j,j),hr);
			for (i=0L;i<=n;i++) erg+=mult_apply(hr,S_M_IJ(b,i,j));
			erg+=copy(hr,S_M_IJ(b,j,j));
			erg+=addinvers_apply(hr);
			for (k=0L;k<=n;k++)
			{
				if (k==j) k++; /* spalte j nicht anwenden */
				if (k>n) break;
				for (i=0L;i<=n;i++)
				{
				if (i==j) i++; /* auf zeile j nicht anwenden */
				if (i>n) break;
				erg+=mult(S_M_IJ(b,i,j),S_M_IJ(b,j,k),hs);

				erg += addinvers_apply(hs);
				erg += add_apply(hs,S_M_IJ(b,i,k));

				};
				erg+=mult_apply(hr,S_M_IJ(b,j,k));
			};
		}; /* end else */
	}; /* end while */
	if (erg != OK)
		goto im_ende;
	erg+=SYMfreeall(hr); 
	erg+=SYMfreeall(hs);

	erg+=m_il_v(n+1L,hv);
	if (not singulaer)
		/*spaltentausch*/
		for (i=0L;i<=n;i++)
			{
			for (k=0L;k<=n;k++) {
				erg+=copy(S_M_IJ(b,i,k), S_V_I(hv,S_V_II(p,k)));
				}
			for (k=0L;k<=n;k++) {
				erg+=copy(S_V_I(hv,k), S_M_IJ(b,i,k));
				}
			}

	erg+=SYMfreeall(p);
	erg+=SYMfreeall(hv);
	if (singulaer)	{ 
		freeself(b);
		error("invers_matrix: singulary");
		return(SINGULAER); 
	};
im_ende:
	ENDR("invers_matrix");

}
#endif /* MATRIXTRUE */


#ifdef MATRIXTRUE
INT transpose_matrix(a,b) OP a,b;
/* AK 280388 */ /* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 210891 V1.3 */
{
	INT i,j;
	INT erg = OK;

	erg += m_ilih_m(S_M_HI(a),S_M_LI(a),b);
	C_O_K(b,S_O_K(a));

	for (i=0;i<S_M_HI(b);i++)
		for (j=0;j<S_M_LI(b);j++)
			erg += copy(S_M_IJ(a,j,i),S_M_IJ(b,i,j));

	ENDR("transpose_matrix");
}
#endif/* MATRIXTRUE */


#ifdef MATRIXTRUE
INT comp_kranztafel(a,b) OP a, b;
/* AK 280390 V1.1 */ /* AK 210891 V1.3 */
{
	INT i,j,erg;
	OP x,y;
	
	x = S_M_S(a);
	y = S_M_S(b);	
	for (i=0L;i<S_M_HI(a);i++)
		{
		if 	(i >= S_M_HI(b)) return(1L);
		else    {
			for (j=(INT)0;j<S_M_LI(a);j++)
				{
				if (j>=S_M_LI(b)) return(1L);
				else {
					erg=COMP_INTEGER_INTEGER(x,y);
					if (erg != (INT)0) return(erg);
					x++;y++;
					};
				}
			}
		}
	return((INT)0); /* matrizen sind gleich */
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT comp_matrix(a,b) OP a, b;
/* AK 060789 V1.0 */ /* AK 070290 V1.1 */ /* AK 210891 V1.3 */
{
	INT i,j,erg;
	OP x,y;
	
	x = S_M_S(a);
	y = S_M_S(b);	
	for (i=(INT)0;i<S_M_HI(a);i++)
		{
		if 	(i >= S_M_HI(b)) return(1L);
		else    {
			for (j=(INT)0;j<S_M_LI(a);j++)
				{
				if (j>=S_M_LI(b)) return(1L);
				else {
					erg=comp(x++,y++);
					if (erg != (INT)0) return(erg);
					};
				}
			}
		}
	if ( S_M_HI(b) > S_M_HI(a) ) return(-1L); /* AK 170790 */
	if ( S_M_LI(b) > S_M_LI(a) ) return(-1L); /* AK 241195 */
	return((INT)0); /* matrizen sind gleich */
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT add_apply_matrix_matrix(a,b) OP a,b;
/* AK 220390 V1.1 */ /* AK 090891 V1.3 */
/* AK 210891 V1.3 */
{
	OP c,d;
	INT i,erg = OK;
	if (
		(S_M_HI(a) == S_M_HI(b))&&
		(S_M_LI(a) == S_M_LI(b))
		)
		{
		i = S_M_HI(a)*S_M_LI(a);
		c = S_M_S(a);d=S_M_S(b);
		while (i-- > 0)
			erg += add_apply(c++,d++);
		}
	else    {
		c = callocobject();
		*c = *b;
		C_O_K(b,EMPTY);
		erg += add_matrix_matrix(a,c,b);
		erg += SYMfreeall(c);
		}
	return erg;
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT add_apply_matrix(a,b) OP a,b;
/* AK 220390 V1.1 */
/* AK 210891 V1.3 */
{
	if (EMPTYP(b)) 
		return(copy_matrix(a,b));
	switch (S_O_K(b)){
		case KRANZTYPUS:
		case MATRIX: return(add_apply_matrix_matrix(a,b));
		default:
			printobjectkind(b);
			error("add_apply_matrix:wrong second type");
			return(ERROR);
	}
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT add_matrix(a,b,c) OP a,b,c;
{
	INT erg=OK;
	if (not MATRIXP(a))
		{ erg += WTO("add_matrix",a);goto endr_ende; }
	if (MATRIXP(b))
		erg += add_matrix_matrix(a,b,c);
	else
		erg += WTO("add_matrix",b);
	ENDR("add_matrix");
}
INT add_matrix_matrix(a,b,ergeb) OP a,b,ergeb;
/* AK 041186 */ /* AK 171186 */ /* AK 060789 V1.0 */ /* AK 081289 V1.1 */
/* AK 090891 V1.3 */
{
	INT i,j;
	OP 	len = callocobject(),
	    height = callocobject();
	OP z;
	INT erg = OK;

	if 	(S_M_LI(a) >=S_M_LI(b))
		COPY_INTEGER(S_M_L(a),len);
	else	COPY_INTEGER(S_M_L(b),len);
	if 	(S_M_HI(a) >=S_M_HI(b))
		COPY_INTEGER(S_M_H(a),height);
	else	COPY_INTEGER(S_M_H(b),height);

	erg += b_lh_m(len,height,ergeb);
	C_O_K(ergeb,S_O_K(a));

	z = S_M_S(ergeb);
	for (i=0;i<S_M_HI(ergeb);i++)
		for (j=0;j<S_M_LI(ergeb);j++,z++)
		{
			if (
			    (i<S_M_HI(a))&&(i<S_M_HI(b))&&
			    (j<S_M_LI(a))&&(j<S_M_LI(b))
			    )

			{ 
			erg += add(S_M_IJ(a,i,j), S_M_IJ(b,i,j), z); 
			}

			else if ( (i<S_M_HI(a))&&(j<S_M_LI(a)))
				erg += copy(S_M_IJ(a,i,j),z); 
			else if ( (i<S_M_HI(b))&&(j<S_M_LI(b)))
				erg += copy(S_M_IJ(b,i,j),z); 
		}
	ENDR("add_matrix_matrix");
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT copy_kranztypus(a,b) OP a,b;
/* AK 270390 V1.1 */ /* AK 210891 V1.3 */
/* AK 040392 cast to char * */
	{
	m_ilih_m(S_M_LI(a),S_M_HI(a),b);
	C_O_K(b,S_O_K(a));
	memcpy((char *) S_M_S(b),(char *) S_M_S(a),
			S_M_LI(a)*S_M_HI(a)*sizeof(struct object));
	return(OK);
	}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT copy_matrix(von,b) OP von , b;
/* AK 051186 */ /* AK 021286 */ /* AK 060789 V1.0 */ /* AK 071289 V1.1 */
/* AK 210891 V1.3 */
{
	INT k;
	OP z,w;
	INT erg = OK;

	erg += m_ilih_m(S_M_LI(von),S_M_HI(von),b);
	C_O_K(b,S_O_K(von));


	z = S_M_IJ(von,S_M_HI(von)-1L,S_M_LI(von)-1L);
	w = S_M_IJ(b,S_M_HI(von)-1L,S_M_LI(von)-1L);
	k = S_M_HI(von) * S_M_LI(von);
	for (;k>(INT)0;k--,z--,w--) 
			if (S_O_K(z) == INTEGER) 
				*w = *z;
			else if (EMPTYP(z)) 
				*w = *z;
			else 
				erg += copy(z,w);

	ENDR("copy_matrix");
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT freeself_kranztypus(a) OP a;
/* AK 270390 V1.1 */ /* AK 210891 V1.3 */
	{
	OBJECTSELF d;
	d=S_O_S(a);
	SYM_free(S_M_S(a)); 
	SYMfreeall(S_M_L(a)); 
	SYMfreeall(S_M_H(a));
	SYM_free(d.ob_matrix);
	C_O_K(a,EMPTY);
	return(OK);
	}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT freeself_matrix(matrix) OP matrix;
/* AK 060789 V1.0 */ /* AK 071289 V1.1 */ /* AK 100691 V1.2 */ 
/* AK 160891 V1.3 */
{
	INT k;
	OBJECTSELF d;
	OP z;
	INT erg = OK;

	d=S_O_S(matrix);

	z = S_M_IJ(matrix,S_M_HI(matrix)-1L,S_M_LI(matrix)-1L);
	k = S_M_HI(matrix) * S_M_LI(matrix);
	for (;k>(INT)0;k--,z--) 
			if (S_O_K(z) == INTEGER) ;
			else if (EMPTYP(z));
			else 
				erg += freeself(z);
	
	SYM_free(S_M_S(matrix)); 
	erg += SYMfreeall(S_M_L(matrix)); 
	erg += SYMfreeall(S_M_H(matrix));
	SYM_free(d.ob_matrix);
	C_O_K(matrix,EMPTY);
	ENDR("freeself_matrix");
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
static struct matrix * callocmatrix()
/* AK 060789 V1.0 */ /* AK 220390 V1.1 */ /* AK 160891 V1.3 */
{
	struct matrix *ergebnis;

	ergebnis = (struct matrix *) SYM_calloc((int)1,sizeof(struct matrix));
	if (ergebnis == NULL)
		no_memory();
	return(ergebnis);
}
#endif/* MATRIXTRUE */

#ifdef MATRIXTRUE
INT scan_integermatrix(ergebnis) OP ergebnis;
/* AK 141293 */
{
	return scan_matrix_co(ergebnis, INTEGER);
}
INT scan_matrix(ergebnis) OP ergebnis;
/* AK 141293 */
{
	return scan_matrix_co(ergebnis, EMPTY);
}


INT scan_skewsymmetric_matrix(ergebnis) OP ergebnis;
/* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 080891 V1.3 */
{
	OP len, height;
	INT i,j;
	char a[20];  /* AK 080891 */
	OBJECTKIND kind;

	len = callocobject();
	height = callocobject();
aaa:
	printeingabe("height of skew symmetric matrix"); 
	scan(INTEGER,height);
	copy(height,len);
	printeingabe("enter kind of matrix elements");
	kind=scanobjectkind();
	if (S_I_I(len) <= (INT)0) /* AK 170795 */
		{
		printeingabe("you entered wrong length (<=0), do it again");
		goto aaa;
		}
	if (S_I_I(height) <= (INT)0) /* AK 170795 */
		{
		printeingabe("you entered wrong height (<=0), do it again");
		goto aaa;
		}
	b_lh_m(len,height,ergebnis);
	for (i=0; i<S_I_I(height); i++)
		m_i_i(0L,S_M_IJ(ergebnis,i,i));
	for (i=0; i<S_I_I(height); i++)
	{
		sprintf(a,"row nr %d \n",(i+1L));  /* AK 080891 */
		printeingabe(a);  /* AK 080891 */
		for (j=i+1;j<S_I_I(len);j++)
			{
			scan(kind,S_M_IJ(ergebnis,i,j));
			addinvers(S_M_IJ(ergebnis,i,j), S_M_IJ(ergebnis,j,i));
			}
	};
	return(OK);
}
static INT scan_matrix_co(ergebnis,kind) OP ergebnis; OBJECTKIND kind;
/* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 080891 V1.3 */
{
	OP len, height;
	INT i,j;
	char a[20];  /* AK 080891 */

	len = callocobject();
	height = callocobject();
aaa:
	printeingabe("height of matrix"); 
	scan(INTEGER,height);
	printeingabe("length of matrix"); 
	scan(INTEGER,len);
	if (kind == EMPTY)
		{
		printeingabe("enter kind of matrix elements");
		kind=scanobjectkind();
		}
	if (S_I_I(len) <= (INT)0) /* AK 170795 */
		{
		printeingabe("you entered wrong length (<=0), do it again");
		goto aaa;
		}
	if (S_I_I(height) <= (INT)0) /* AK 170795 */
		{
		printeingabe("you entered wrong height (<=0), do it again");
		goto aaa;
		}
	b_lh_m(len,height,ergebnis);
	for (i=0; i<S_I_I(height); i++)
	{
		sprintf(a,"row nr %d \n",(i+1L));  /* AK 080891 */
		printeingabe(a);  /* AK 080891 */
		for (j=0;j<S_I_I(len);j++)
			scan(kind,S_M_IJ(ergebnis,i,j));
	};
	return(OK);
}
#endif/* MATRIXTRUE */

#ifdef MATRIXTRUE
INT change_column_ij(a,i,j) OP a; INT i,j;
/* AK 301288 vertauscht spalten i und j dabei kein copy */
/* AK 060789 V1.0 */ /* AK 050390 V1.1 */ /* AK 160891 V1.3 */
{
	INT k,erg = OK;
	if (not matrixp(a))
		{
		erg += WTO("change_column_ij",a);
		goto a1;
		}
	if ((i<0) || (i>=S_M_LI(a))) 
		{
		erg += error("change_column_ij: wrong i");
		goto a1;
		}
	if ((j<0) || (j>=S_M_LI(a))) 
		{
		erg += error("change_column_ij: wrong j");
		goto a1;
		}

	for (k=0; k<S_M_HI(a); k++) 
		erg += swap(S_M_IJ(a,k,i),S_M_IJ(a,k,j));
a1:
	ENDR("change_column_ij");
}
#endif/* MATRIXTRUE */



#ifdef MATRIXTRUE
INT change_row_ij(a,i,j) OP a; INT i,j;
/* AK 301288 vertauscht zeilen i und j dabei kein copy */
/* AK 060789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{
	INT k,erg = OK;

	if (not matrixp(a))
		{
		erg += WTO("change_row_ij",a);
		goto a1;
		}
	if ((i<(INT)0) || (i>=S_M_HI(a))) 
		{
		fprintf(stderr,"i = %ld, S_M_HI = %ld\n",i, S_M_HI(a));
		erg += error("change_row_ij: wrong i");
		goto a1;
		}
	if ((j<(INT)0) || (j>=S_M_HI(a))) 
		{
		fprintf(stderr,"j = %ld, S_M_HI = %ld\n",j, S_M_HI(a));
		erg += error("change_row_ij: wrong j");
		goto a1;
		}

	for (k=(INT)0; k<S_M_LI(a); k++) 
		erg += swap(S_M_IJ(a,i,k),S_M_IJ(a,j,k));
a1:
	if (erg != OK)
		EDC("change_row_ij");
	return erg;
}
#endif/* MATRIXTRUE */

#ifdef MATRIXTRUE
OP s_m_s(a) OP a;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ OBJECTSELF c; c = s_o_s(a); return(c.ob_matrix->m_self); }

OP s_m_h(a) OP a;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ OBJECTSELF c; c = s_o_s(a); return(c.ob_matrix->m_height); }

OP s_m_l(a) OP a;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ OBJECTSELF c; c = s_o_s(a); return(c.ob_matrix->m_length); }

INT s_m_hi(a) OP a;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ return(s_i_i(s_m_h(a))); }

INT s_m_li(a) OP a;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ return(s_i_i(s_m_l(a))); }

INT c_m_s(a,b) OP a,b;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ OBJECTSELF c; c = s_o_s(a); c.ob_matrix->m_self = b; return(OK); }


INT c_m_h(a,b) OP a,b;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ OBJECTSELF c; c = s_o_s(a); c.ob_matrix->m_height = b; return(OK); }


INT c_m_l(a,b) OP a,b;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ OBJECTSELF c; c = s_o_s(a); c.ob_matrix->m_length = b; return(OK); }

OP s_m_ij(a,i,j) OP a; INT i,j;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ 
	if (not MATRIXP(a))
		{
		printobjectkind(a);
		error("s_m_ij:no matrix object");
		}
	if (i < (INT)0)
		{
		debugprint(a);
		fprintf(stderr,"index = %ld\n",i);
		error("s_m_ij:row index too small");
		}
	if (i >= s_m_hi(a))
		{
		debugprint(a);
		fprintf(stderr,"index = %ld\n",i);
		error("s_m_ij:row index too big");
		}
	if (j >= s_m_li(a))
		{
		debugprint(a);
		fprintf(stderr,"index = %ld\n",j);
		error("s_m_ij:column index too big");
		}
	if (j < (INT)0)
		{
		debugprint(a);
		fprintf(stderr,"index = %ld\n",j);
		error("s_m_ij:column index too small");
		}
	return(s_m_s(a) + (s_m_li(a)*i+j) ); 
}

INT s_m_iji(a,i,j) OP a; INT i,j;
/* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
{ return(s_i_i(s_m_ij(a,i,j))); }
#endif /* MATRIXTRUE */


#ifdef MATRIXTRUE
INT fprint_matrix(f,obj) FILE  *f; OP obj;
/* AK 211186 */ /* AK 070789 V1.0 */ /* AK 181289 V1.1 */
/* AK 210891 V1.3 */
{
	INT i,j;
	for (i=0;i<S_M_HI(obj);i++)
	{
		fprintf(f,"\n[");
		if (f == stdout) zeilenposition=0;
		for (j=0;j<S_M_LI(obj);j++)
		{
			fprint(f,S_M_IJ(obj,i,j));
			fprintf(f,":");zeilenposition++;
			if ((f == stdout)&&(zeilenposition>70L))
				{fprintf(stdout,"\n");zeilenposition = 0L;}
		};
		fprintf(f,"]");
	};
	fprintf(f,"\n");
	if (f == stdout) zeilenposition=0L;
	return(OK);
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT tex_matrix(obj) OP obj;
{
	return tex_matrix_co(obj,tex);
}

INT tex_matrix_co(obj,f) OP obj; INT (*f)();
/* AK 150988 */ /* AK 310790 V1.1 */
/* AK 070291 V1.2 texout for output */ /* AK 210891 V1.3 */
{
	INT i,j;
	INT ts = texmath_yn; /* AK 190892 */
	INT erg = OK;

	fprintf(texout,"\n");
	if (texmath_yn == 0L)  /* d.h. not in math mode */
		{
		fprintf(texout,"$");
		texmath_yn=1L;
		}
	fprintf(texout,"\\matrix { \n");
	texposition = 0L;
	for (i=0;i<S_M_HI(obj);i++)
	{
		for (j=0L;j<S_M_LI(obj);j++)
		{
			(*f)(S_M_IJ(obj,i,j));
			fprintf(texout," & "); 
			texposition += 3L;
		}
		fprintf(texout," \\cr\n");
		texposition=0L;
	};
	fprintf(texout," }");
	if (ts == 0L) 
		{
		fprintf(texout,"$");
		texmath_yn=0L;
		}
	fprintf(texout," \n");
	texposition=0L;
	ENDR("tex_matrix_co");
}


INT mult_scalar_matrix(scalar,mat,res) OP scalar , mat,res;
/* AK 310588 */ /* komponentenweise multiplikation */
/* AK 010889 V1.0 */ /* AK 081289 V1.1 */ /* AK 210891 V1.3 */
{
	INT i,j,erg = OK;
	OP height,len;

	if (S_M_S(mat) == NULL) /* AK 161092 */
		return error("mult_scalar_matrix:self == NULL");
	height = callocobject();
	len = callocobject();

	COPY_INTEGER(S_M_L(mat), len); 
	COPY_INTEGER(S_M_H(mat), height);
	erg += b_lh_m(len,height,res);
	for (i=0L;i<S_I_I(height);i++) 
		for(j=0L;j<S_I_I(len);j++)
			erg += mult(scalar,S_M_IJ(mat,i,j),S_M_IJ(res,i,j));
	return erg;
}
#endif/* MATRIXTRUE */

#ifdef MATRIXTRUE
INT mult_apply_scalar_matrix(a,b) OP a,b;
/* AK 150290 V1.1 */ /* AK 210891 V1.3 */
{
	OP z = S_M_S(b);
	INT grenze = S_M_LI(b)*S_M_HI(b);
	INT i;
	for (i=0L;i<grenze;i++,z++) mult_apply(a,z);
	return(OK);
}
#endif/* MATRIXTRUE */

#ifdef MATRIXTRUE
INT mult_apply_matrix_matrix(a,b) OP a,b;
/* AK 131190 V1.1 */ /* b =  a * b */
/* AK 210891 V1.3 */
{
	OP c = callocobject();
	INT erg = OK; /* AK 200192 */
	*c = *b;
	C_O_K(b,EMPTY);
	erg += mult_matrix_matrix(a,c,b);
	erg += SYMfreeall(c);
	ENDR("mult_apply_matrix_matrix");
}
#endif/* MATRIXTRUE */

#ifdef MATRIXTRUE
INT mult_matrix_matrix(a,b,c) OP a,b,c;
/* AK 280388 */ /* AK 060789 V1.0 */ /* AK 111289 V1.1 */
/* c = a * b */
/* AK 210891 V1.3 */
{
	INT i,j,k;
	OP z; /* zwischen ergebnis bei matrix-multiplikation */

	if (neq(s_m_l(a),s_m_h(b)))
	{ 
		error("mult_matrix_matrix:can not be multiplied");
		return(ERROR); 
	};
	m_ilih_m(S_M_LI(b),S_M_HI(a),c);
	z=callocobject(); /* zwischensumme*/
	for (i=0;i<S_M_HI(a);i++)	/* ueber zeilen der linken Matrix */
		for (j=0;j<S_M_LI(b);j++) /* ueber spalten der rechten Matrix */
			for (k=0;k<S_M_LI(a);k++)
			{ 
				mult(S_M_IJ(a,i,k),S_M_IJ(b,k,j),z);
				add(z,S_M_IJ(c,i,j),S_M_IJ(c,i,j)); 
			};
	SYMfreeall(z); 
	return(OK);
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT mult_matrix(a,b,d) OP a,b,d;
/* AK 070789 V1.0 */ /* AK 310790 V1.1 */ /* AK 210891 V1.3 */
{
	switch(S_O_K(b))
	{
	case BRUCH:
	case INTEGER:
	case LONGINT: return mult_scalar_matrix(b,a,d);
	case MATRIX: return mult_matrix_matrix(a,b,d);
	case VECTOR: return mult_matrix_vector(a,b,d); /* AK 200192 */
	default:
		{
			printobjectkind(b);
			error("mult_matrix:wrong second type"); 
			return(ERROR);
		}
	}
}
#endif/* MATRIXTRUE */

#ifdef MATRIXTRUE
#ifdef VECTORTRUE
INT mult_matrix_vector(b,a,c) OP a, b, c;
/* AK 200192 */
{
	INT i,j;
	INT erg = OK;
	OP d;
	if (S_O_K(a) != VECTOR) return ERROR;
	if (S_O_K(b) != MATRIX) return ERROR;
	if ((a == c) || (b == c)) return ERROR;
	if (S_V_LI(a) != S_M_LI(b)) return ERROR;
	erg += m_il_nv(S_M_HI(b),c);
	d = callocobject();
	for (i=0L;i<S_V_LI(c);i++)
	for (j=0L;j<S_V_LI(a);j++)
		{
		erg += mult(S_M_IJ(b,i,j),S_V_I(a,j),d);
		erg += add_apply(d,S_V_I(c,i));
		}
	erg += SYMfreeall(d);
	ENDR("mult_matrix_vector");
}
#endif /* VECTORTRUE */
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT mult_apply_matrix(a,b) OP a,b;
/* AK 131190 V1.1 */ /* AK 210891 V1.3 */
{
	switch(S_O_K(b))
	{
	case MATRIX: return(mult_apply_matrix_matrix(a,b));
	default:
		{
			printobjectkind(b);
			error("mult_apply_matrix:wrong second type"); 
			return(ERROR);
		}
	}
}
#endif/* MATRIXTRUE */


#ifdef MATRIXTRUE
INT objectread_matrix(fp,matrix) FILE *fp; OP matrix;
/* AK 300888 */ /* AK 070789 V1.0 */ /* AK 310790 V1.1 */
/* AK 210891 V1.3 */
{
	INT i,j;
	OP l= callocobject();
	OP h = callocobject();
	objectread(fp,h);
	objectread(fp,l);
	b_lh_m(l,h,matrix);
	for (i=0;i<S_M_HI(matrix); i++)
		for (j=0;j<S_M_LI(matrix); j++)
			objectread(fp,S_M_IJ(matrix,i,j));
	return(OK);
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT objectwrite_matrix(fp,matrix) FILE *fp; OP matrix;
/* AK 300888 */ /* AK 070789 V1.0 */ /* AK 310790 V1.1 */
/* AK 210891 V1.3 */
{
	INT i,j;

	fprintf(fp, " %d ",MATRIX);
	objectwrite(fp,S_M_H(matrix));
	/* zuerst die hoehe */
	objectwrite(fp,S_M_L(matrix));
	/* dann die laenge */

	for (i=0;i<S_M_HI(matrix); i++)
		for (j=0;j<S_M_LI(matrix); j++)
			objectwrite(fp,S_M_IJ(matrix,i,j));
	return(OK);
}
#endif/* MATRIXTRUE */

#ifdef MATRIXTRUE
INT test_matrix() 
/* AK 181289 V1.1 */
/* AK 120891 V1.3 */
{
	OP a = callocobject();
	OP b = callocobject();

	printf("test_matrix:scan(a)");
	scan(MATRIX,a);
	println(a);
	printf("test_matrix:add(a,a,b)");
	add(a,a,b);
	println(b);
	printf("test_matrix:mult(a,b,b)"); 
	mult(a,b,b); 
	println(b);
	printf("test_matrix:kronecker_product(a,b,b)"); 
	kronecker_product(a,b,b); 
	println(b);
#ifdef BRUCHTRUE
	printf("test_matrix:invers(b,a)"); 
	invers(b,a); 
	println(a);
#endif /* BRUCHTRUE */
	printf("test_matrix:delete_row_matrix(a,1L,b)");
	delete_row_matrix(a,1L,b); 
	println(b);
	printf("test_matrix:delete_column_matrix(b,1L,b)");
	delete_column_matrix(b,1L,b); 
	println(b);


	SYMfreeall(a);
	SYMfreeall(b);
	return(OK);
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT trace_matrix(a,b) OP a,b;
/* AK 131289 spur einer matrix */
/* AK 131289 V1.1 */ /* AK 270291 V1.2 */ /* AK 210891 V1.3 */
{
	INT i,erg=OK;
	if (not MATRIXP(a)) 
		return WTO(a,"trace_matrix");
	if (not quadraticp(a)) 
		return error("trace_matrix: matrix not quadratic");
	erg += m_i_i(0L,b);
	for (i=S_M_HI(a)-1L; i>=0L; i--)
		erg += add_apply(S_M_IJ(a,i,i),b);
	if (erg != OK)
		EDC("trace_matrix");
	return erg;
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT symmetricp_matrix(a) OP a;
/* AK 150296 */
{
	INT i,j;
	if (S_M_HI(a) != S_M_LI(a)) return FALSE;
	for (i=0L;i<S_M_HI(a) ; i++)
	for (j=0L;j<i;j++)
		if (neq(S_M_IJ(a,i,j),S_M_IJ(a,j,i))) return FALSE;
	return TRUE;
}

#ifdef VECTORTRUE
INT spalten_summe(a,b) OP a,b;
/* AK 281289 summe ueber die Spalten, ergebnis ist vector */
/* AK 281289 V1.1 */ /* AK 240791 V1.3 */
{
	INT i,j;
	INT erg = OK;
	if (not MATRIXP(a))
		return(error("spalten_summe: typ not MATRIX"));
	erg += m_il_nv(S_M_LI(a),b);
	for (j=0;j<S_M_LI(a);j++)
		for (i=0;i<S_M_HI(a);i++)
			erg += add_apply(S_M_IJ(a,i,j),S_V_I(b,j));
	ENDR("spalten_summe");
}
#endif /* VECTORTRUE */
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
char * t_INTMATRIX_charvektor(a) OP a;
/* AK 210891 V1.3 */
{
	INT i,j,k=0;
	char *erg = (char *) 
		SYM_malloc(S_M_HI(a)*S_M_LI(a)*sizeof(char));
	if(erg == NULL) error("t_INTMATRIX_charvektor:no memory");
	else {
		for (i=0;i<S_M_HI(a);i++)
				for (j=0;j<S_M_LI(a);j++,k++)
					erg[k] =(char)S_M_IJI(a,i,j);
		}
	return(erg);
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
#ifdef VECTORTRUE
INT m_vector_diagonalmatrix(a,b) OP a,b;
/* AK 171290 V1.1 */ /* AK 110691 V1.2 */ /* AK 210891 V1.3 */
{
	INT i;
	m_lh_nm(S_V_L(a),S_V_L(a),b);
	for (i=0L; i<S_M_HI(b);i++)
			copy(S_V_I(a,i),S_M_IJ(b,i,i));
	return OK;
}
#endif /* VECTORTRUE */
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT max_matrix(a,b) OP a,b;
/* b becomes copy of the maximum entry */
/* AK 110691 V1.2 */ /* AK 210891 V1.3 */
{
	OP z = S_M_S(a),zb = S_M_S(a);
	INT i;

	for (i=S_M_HI(a)*S_M_LI(a); i>0L; i--, z++)
		if (GR(z,zb)) zb = z;
	return copy(zb,b);
}

INT zeilen_summe(a,b) OP a,b;
/* AK 281289 summe ueber die Zeilen, ergebnis ist vector */
/* AK 281289 V1.1 */ /* AK 210891 V1.3 */
{
	INT i,j;
	INT erg = OK;
	if (S_O_K(a) != MATRIX) 
		return(error("spalten_summe: typ not MATRIX"));
	
	erg += m_il_nv(S_M_HI(a),b);
	for (j=0;j<S_M_HI(a);j++)
		for (i=0;i<S_M_LI(a);i++)
			erg += add_apply(S_M_IJ(a,j,i),S_V_I(b,j));
	ENDR("zeilen_summe");
}


INT n_fold_kronecker_product(n,a,b) OP	n; OP	a; OP	b;
/* RH 080891 V1.3 */
{
	INT	i;

	if(S_I_I(n)>= 2L) 	kronecker_product(a,a,b);
	for(i=2;i<S_I_I(n);++i) kronecker_product(a,b,b);
	return OK;
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE

INT kronecker_product(a,b,c) OP a,b,c;
/* RH 080891 V1.3 */
{
	INT	i;
	INT	j;
	INT	k;
	INT	l,erg=OK;

	OP	H,dim;

        if (check_equal_3(a,b,c,kronecker_product,&erg) == EQUAL)
                goto ke;


	H 	=	callocobject();
	dim	=	callocobject();

	erg += mult(S_M_L(a),S_M_L(b),dim);

	erg += m_lh_m(dim,dim,c);
	erg += m_lh_m(S_M_L(a),S_M_L(b),H);

	for(i=0L;i<S_M_LI(a);++i)
	{
		for(j=0L;j<S_M_HI(a);++j)
		{
			erg += mult(S_M_IJ(a,i,j),b,H);
			for(k=0L;k<S_M_LI(H);++k)
				for(l=0L;l<S_M_HI(H);++l)
			erg += copy(S_M_IJ(H,k,l),
			 S_M_IJ(c,i*S_M_HI(b)+k,j*S_M_LI(H)+l));
		}
	}

	SYMfreeall(H);
	SYMfreeall(dim);
ke:
	return erg;
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE

INT select_row(a,index,b) OP a,b; INT index;
/* transfer one row to vector object  */
/* AK 110592 */
{
	INT j;
	INT erg = OK;
	if (tableauxp(a)) /* AK 280193 */
		return select_row_tableaux(a,index,b);
	if (not matrixp(a))
		return error("select_row:only for MATRIX");
	erg += m_il_v(S_M_LI(a),b);
	for (j=0L;j<S_M_LI(a);j++)
		erg += copy(S_M_IJ(a,index,j),S_V_I(b,j));
	return erg;
}
#endif /* MATRIXTRUE */

#ifdef MATRIXTRUE
INT select_column(a,index,b) OP a,b; INT index;
/* transfer one column to vector object  */
/* AK 110592 */
{
	INT j;
	INT erg = OK;
	if (tableauxp(a)) /* AK 280193 */
		return select_column_tableaux(a,index,b);
	if (not matrixp(a))
		return error("select_column:only for MATRIX");
	erg += m_il_v(S_M_HI(a),b);
	for (j=0L;j<S_M_HI(a);j++)
		erg += copy(S_M_IJ(a,j,index),S_V_I(b,j));
	return erg;
}


INT rank(a,b) OP a,b;
/* AK 220294 */
{
	INT i,j,k,l,z=0;
	INT res = 0, erg = OK;
	OP e,c,d;

    	if (check_equal_2(a,b,rank,&erg) == EQUAL)
               return erg;


	c = callocobject();
	d = callocobject();
	e = callocobject();
	erg += copy(a,e);
	a = e;
	j=(INT)0;
a2:
	for (i=(INT)0;i<S_M_HI(a);i++)
		{
		if (not nullp(S_M_IJ(a,i,j))) 
			goto a1;
		}
	j++;
	if (j == S_M_LI(a))
		goto a3;
	goto a2;
a1:
	/* i , j ist nicht null */
	if (i != z) 
		erg += change_row_ij(a,i,z); /* AK 080394 */
	i = z++;

	res ++;
	for (k = (INT)0;k<S_M_HI(a);k++)
		if (k != i)
			{
			erg += div3(S_M_IJ(a,k,j),S_M_IJ(a,i,j),c);
			for (l=(INT)0;l<S_M_LI(a);l++)
				{
				erg += mult(c,S_M_IJ(a,i,l),d);
				erg += addinvers_apply(d);
				erg += add_apply(d,S_M_IJ(a,k,l));
	/*
				erg += sub(S_M_IJ(a,k,l),d,S_M_IJ(a,k,l));
	*/
				}
			}
	for (l=j+1;l<S_M_LI(a);l++)
		/*
		erg += sub(S_M_IJ(a,i,l),S_M_IJ(a,i,l),S_M_IJ(a,i,l));
		*/
		m_i_i((INT)0,S_M_IJ(a,i,l));
	j++;
	if (j != S_M_LI(a))
		goto a2;
a3:
	erg += m_i_i(res,b);
	erg += SYMfreeall(c);	
	erg += SYMfreeall(d);	
	erg += SYMfreeall(a);	
	ENDR("rank");
}

#endif /* MATRIXTRUE */




