#include <stdio.h>
#include <string.h>
#include <math.h>
#include <limits.h>

#include "output.h"
#include "funcs.h"
#include "matrix.h"
#include "linear.h"
#include "interval.h"
#include "mainloop.h"

#define wrong_arg() { error=26; output("Wrong argument\n"); return; }

void minmax (double *x, long n, double *min, double *max, 
	int *imin, int *imax)
/***** minmax
	compute the total minimum and maximum of n double numbers.
*****/
{	long i;
	if (n==0)
	{	*min=0; *max=0; *imin=0; *imax=0; return; }
	*min=*x; *max=*x; *imin=0; *imax=0; x++;
	for (i=1; i<n; i++)
	{	if (*x<*min) { *min=*x; *imin=(int)i; }
		else if (*x>*max) { *max=*x; *imax=(int)i; }
		x++;
	}
}

/****************************************************************************
 *	utility functions
 *
 ****************************************************************************/
void msize (header *hd)
{	header *result,*st=hd,*hd1,*end=(header *)newram;
	int r,c,r0=0,c0=0;
	if (!hd) wrong_arg_in("size");
	result=new_matrix(1,2,""); if (error) return;
	while (end>hd)
	{	hd1=getvariable(hd); if (!hd1) varnotfound("size");
		if (hd1->type==s_matrix || hd1->type==s_cmatrix
			|| hd1->type==s_imatrix)
		{	r=dimsof(hd1)->r;
			c=dimsof(hd1)->c;
		}
		else if (hd1->type==s_real || hd1->type==s_complex
			|| hd1->type==s_interval)
		{	r=c=1;
		}
		else if (hd1->type==s_submatrix || hd1->type==s_csubmatrix
			|| hd1->type==s_isubmatrix)
		{	r=submdimsof(hd1)->r;
			c=submdimsof(hd1)->c;
		}
		else wrong_arg_in("size");
		if ((r>1 && r0>1 && r!=r0) || (c>1 && c0>1 && c!=c0))
		{	if (r0!=r && c0!=c)
			{	output("Matrix dimensions must agree for size!\n");
				error=1021; return;
			}
		}
		else
		{	if (r>r0) r0=r;
			if (c>c0) c0=c;
		}
		hd=nextof(hd);
	}
	*matrixof(result)=r0;
	*(matrixof(result)+1)=c0;
	moveresult(st,result);
}

void mcols (header *hd)
{	header *st=hd,*res;
	int n;
	hd=getvalue(hd); if (error) return;
	switch (hd->type)
	{	case s_matrix :
		case s_cmatrix :
		case s_imatrix : n=dimsof(hd)->c; break;
		case s_submatrix :
		case s_csubmatrix :
		case s_isubmatrix : n=submdimsof(hd)->c; break;
		case s_real :
		case s_complex :
		case s_interval : n=1; break;
		case s_string : n=(int)strlen(stringof(hd)); break;
		default : wrong_arg_in("cols");
	}
	res=new_real(n,""); if (error) return;
	moveresult(st,res);
}

void mrows (header *hd)
{	header *st=hd,*res;
	int n;
	hd=getvalue(hd); if (error) return;
	switch (hd->type)
	{	case s_matrix :
		case s_cmatrix :
		case s_imatrix : n=dimsof(hd)->r; break;
		case s_submatrix :
		case s_csubmatrix :
		case s_isubmatrix : n=submdimsof(hd)->r; break;
		case s_real :
		case s_complex :
		case s_interval : n=1; break;
		default : wrong_arg_in("rows");
	}
	res=new_real(n,""); if (error) return;
	moveresult(st,res);
}

void mredim (header *hd)
{	header *st=hd,*hd1,*result;
	int c1,r1;
	double *m;
	unsigned long i,n,size1,size;
	hd1=nextof(hd);
	hd=getvalue(hd); if (error) return;
	hd1=getvalue(hd1); if (error) return;
	if (hd1->type!=s_matrix || dimsof(hd1)->r!=1 || dimsof(hd1)->c!=2
		|| (hd->type!=s_matrix && hd->type!=s_cmatrix))
		wrong_arg_in("redim");
	m=matrixof(hd1);
	if (*m<1 || *m>INT_MAX) wrong_arg_in("redim");
	r1=(int)(*m++);
	if (*m<1 || *m>INT_MAX) wrong_arg_in("redim");
	c1=(int)(*m);
	size1=(long)c1*r1;
	size=(long)dimsof(hd)->c*dimsof(hd)->r;
	if (size<size1) n=size;
	else n=size1;
	if (hd->type==s_matrix)
	{	result=new_matrix(r1,c1,""); if (error) return;
		memmove((char *)matrixof(result),(char *)matrixof(hd),
			n*sizeof(double));
		if (n<size1)
		{	m=matrixof(result)+n;
			for (i=n; i<size1; i++) *m++=0.0;
		}
	}
	else if (hd->type==s_cmatrix)
	{	result=new_cmatrix(r1,c1,""); if (error) return;
		memmove((char *)matrixof(result),(char *)matrixof(hd),
			2*n*sizeof(double));
		if (n<size1)
		{	m=matrixof(result)+2*n;
			for (i=n; i<size1; i++) { *m++=0.0; *m++=0.0; }
		}
	}
	else wrong_arg_in("redim");
	moveresult(st,result);
}

void mresize (header *hd)
{	header *st=hd,*hd1,*result;
	int c1,r1,c0,r0,i,j;
	double *m,*m0,*mr,*mm;
	hd1=nextof(hd);
	hd=getvalue(hd); if (error) return;
	hd1=getvalue(hd1); if (error) return;
	if (hd1->type!=s_matrix || dimsof(hd1)->r!=1 || dimsof(hd1)->c!=2)
		wrong_arg_in("redim");
	m=matrixof(hd1);
	if (*m<1 || *m>INT_MAX) wrong_arg_in("redim");
	r1=(int)(*m++);
	if (*m<1 || *m>INT_MAX) wrong_arg_in("redim");
	c1=(int)(*m);
	getmatrix(hd,&r0,&c0,&m0); mm=m0;
	if ((r0!=1 && r1!=r0) || (c0!=1 && c1!=c0))
	{	output("Cannot resize these!\n"); error=1; return;
	}
	if (hd->type==s_matrix || hd->type==s_real)
	{   result=new_matrix(r1,c1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r1; i++)
		{	if (c0==1)
			{	for (j=0; j<c1; j++) *mr++=*m0;
				m0++;
			}
			else for (j=0; j<c1; j++) *mr++=*m0++;
			if (r0==1) m0=mm;
		}
	}
	else if (hd->type==s_cmatrix || hd->type==s_complex)
	{   result=new_cmatrix(r1,c1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r1; i++)
		{	if (c0==1)
			{	for (j=0; j<c1; j++) { *mr++=*m0; *mr++=*(m0+1); }
				m0+=2;
			}
			else for (j=0; j<c1; j++) { *mr++=*m0++; *mr++=*m0++; }
			if (r0==1) m0=mm;
		}
	}
	else if (hd->type==s_imatrix || hd->type==s_interval)
	{   result=new_imatrix(r1,c1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r1; i++)
		{	if (c0==1)
			{	for (j=0; j<c1; j++) { *mr++=*m0; *mr++=*(m0+1); }
				m0+=2;
			}
			else for (j=0; j<c1; j++) { *mr++=*m0++; *mr++=*m0++; }
			if (r0==1) m0=mm;
		}
	}
	else wrong_arg_in("redim");
	moveresult(st,result);
}

/****************************************************************************
 *	generate matrix
 *
 ****************************************************************************/
void vectorize (header *init, header *step, header *end)
{	double vinit,vstep,vend,*m;
	int count;
	header *result,*st=init;
	init=getvalue(init); step=getvalue(step); end=getvalue(end);
	if (error) return;
	if (init->type!=s_real || step->type!=s_real || end->type!=s_real)
	{	output("The : allows only real arguments!\n");
		error=15; return;
	}
	vinit=*realof(init); vstep=*realof(step); vend=*realof(end);
	if (vstep==0)
	{	output("A step size of 0 is not allowed in : \n");
		error=401; return;
	}
	if (fabs(vend-vinit)/fabs(vstep)+1+epsilon>INT_MAX)
	{	output1("A vector can only have %d elements\n",INT_MAX);
		error=402; return;
	}
	count=(int)(floor(fabs(vend-vinit)/fabs(vstep)+1+epsilon));
	if ((vend>vinit && vstep<0) || (vend<vinit && vstep>0))
		count=0;
	result=new_matrix(1,count,""); if (error) return;
	m=matrixof(result);
	while (count>=0)
	{	*m++=vinit;
		vinit+=vstep;
		count--;
	}
	moveresult(st,result);
}

void mmatrix (header *hd)
{	header *st=hd,*hd1,*result;
	long i,n;
	double x,xi;
	double *m,*mr;
	int c,r,c1,r1;
	hd1=nextof(hd);
	hd=getvalue(hd);
	if (error) return;
	hd1=getvalue(hd1);
	if (error) return;
	if (hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		if (*m<0 || *m>INT_MAX || *(m+1)<0 || *(m+1)>INT_MAX)
			wrong_arg_in("matrix");
		r1=(int)*m; c1=(int)*(m+1);
		if (hd1->type==s_real)
		{	result=new_matrix(r1,c1,"");
			mr=matrixof(result);
			x=*realof(hd1);
			n=(long)c1*r1;
			for (i=0; i<n; i++) *mr++=x;
		}
		else if (hd1->type==s_complex)
		{	result=new_cmatrix(r1,c1,"");
			mr=matrixof(result);
			x=*realof(hd1); xi=*(realof(hd1)+1);
			n=(long)c1*r1;
			for (i=0; i<n; i++) 
			{	*mr++=x; *mr++=xi;
			}
		}
		else wrong_arg_in("matrix");
	}
	else wrong_arg_in("matrix");
	moveresult(st,result);
}

void mzerosmat (header *hd)
{	header *result,*st=hd;
	double rows,cols,*m;
	int r,c;
	long i,n;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2)
		wrong_arg_in("zeros");
	rows=*matrixof(hd); cols=*(matrixof(hd)+1);
	if (rows<0 || rows>=INT_MAX || cols<0 || cols>=INT_MAX)
		wrong_arg_in("zeros");
	r=(int)rows; c=(int)cols;
	result=new_matrix(r,c,""); if (error) return;
	m=matrixof(result);
	n=c*r;
	for (i=0; i<n; i++) *m++=0.0;
	moveresult(st,result);
}

void mones (header *hd)
{	header *result,*st=hd;
	double rows,cols,*m;
	int r,c;
	long i,n;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2)
		wrong_arg_in("ones");
	rows=*matrixof(hd); cols=*(matrixof(hd)+1);
	if (rows<0 || rows>=INT_MAX || cols<0 || cols>=INT_MAX)
		wrong_arg_in("ones");
	r=(int)rows; c=(int)cols;
	result=new_matrix(r,c,""); if (error) return;
	m=matrixof(result);
	n=c*r;
	for (i=0; i<n; i++) *m++=1.0;
	moveresult(st,result);
}

void mdiag (header *hd)
{	header *result,*st=hd,*hd1,*hd2=0;
	double rows,cols,*m,*md;
	int r,c,i,ik=0,k,rd,cd;
	long l,n;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_matrix || dimsof(hd)->r!=1 || dimsof(hd)->c!=2)
		wrong_arg_in("diag");
	rows=*matrixof(hd); cols=*(matrixof(hd)+1);
	if (rows<0 || rows>=INT_MAX || cols<0 || cols>=INT_MAX)
		wrong_arg_in("diag");
	r=(int)rows; c=(int)cols;
	hd1=next_param(st); if (hd1) hd2=next_param(hd1);
	if (hd1) hd1=getvalue(hd1);
	if (hd2) hd2=getvalue(hd2);
	if (error) return;
	if	(hd1->type!=s_real) wrong_arg_in("diag");
	k=(int)*realof(hd1);
	if (hd2->type==s_matrix || hd2->type==s_real)
	{	result=new_matrix(r,c,""); if (error) return;
		m=matrixof(result);
		n=(long)c*r;
		for (l=0; l<n; l++) *m++=0.0;
		getmatrix(hd2,&rd,&cd,&md);
		if (rd!=1 || cd<1) wrong_arg_in("diag");
		m=matrixof(result);
		for (i=0; i<r; i++)
		{	if (i+k>=0 && i+k<c)
			{	*mat(m,c,i,i+k)=*md;
				ik++; if (ik<cd) md++;
			}
		}
	}
	else if (hd2->type==s_cmatrix || hd2->type==s_complex)
	{	result=new_cmatrix(r,c,""); if (error) return;
		m=matrixof(result);
		n=(long)2*(long)c*r;
		for (l=0; l<n; l++) *m++=0.0;
		getmatrix(hd2,&rd,&cd,&md);
		if (rd!=1 || cd<1) wrong_arg_in("diag");
		m=matrixof(result);
		for (i=0; i<r; i++)
		{	if (i+k>=0 && i+k<c)
			{	*cmat(m,c,i,i+k)=*md;
				*(cmat(m,c,i,i+k)+1)=*(md+1);
				ik++; if (ik<cd) md+=2;
			}
		}
	}
	else if (hd2->type==s_imatrix || hd2->type==s_interval)
	{	result=new_imatrix(r,c,""); if (error) return;
		m=matrixof(result);
		n=(long)2*(long)c*r;
		for (l=0; l<n; l++) *m++=0.0;
		getmatrix(hd2,&rd,&cd,&md);
		if (rd!=1 || cd<1) wrong_arg_in("diag");
		m=matrixof(result);
		for (i=0; i<r; i++)
		{	if (i+k>=0 && i+k<c)
			{	*imat(m,c,i,i+k)=*md;
				*(imat(m,c,i,i+k)+1)=*(md+1);
				ik++; if (ik<cd) md+=2;
			}
		}
	}
	else wrong_arg_in("diag");
	moveresult(st,result);
}

void msetdiag (header *hd)
{	header *result,*st=hd,*hd1,*hd2=0;
	double *m,*md,*mhd;
	int r,c,i,ik=0,k,rd,cd;
	hd=getvalue(st); if (error) return;
	if (hd->type!=s_matrix && hd->type!=s_cmatrix)
		wrong_arg_in("setdiag");
	getmatrix(hd,&c,&r,&mhd);
	hd1=next_param(st); if (hd1) hd2=next_param(hd1);
	if (hd1) hd1=getvalue(hd1);
	if (hd2) hd2=getvalue(hd2);
	if (error) return;
	if	(hd1->type!=s_real) wrong_arg_in("setdiag");
	k=(int)*realof(hd1);
	if (hd->type==s_matrix &&
			(hd2->type==s_complex || hd2->type==s_cmatrix))
		{	make_complex(st); msetdiag(st); return;
		}
	else if (hd->type==s_cmatrix &&
			(hd2->type==s_real || hd2->type==s_matrix))
		{	make_complex(nextof(nextof(st))); msetdiag(st); return;
		}
	else if (hd->type==s_imatrix &&
			(hd2->type==s_real || hd2->type==s_matrix))
		{	make_interval(nextof(nextof(st))); msetdiag(st); return;
		}
	if (hd->type==s_matrix)
	{	result=new_matrix(r,c,""); if (error) return;
		m=matrixof(result);
		memmove((char *)m,(char *)mhd,(long)c*r*sizeof(double));
		getmatrix(hd2,&rd,&cd,&md);
		if (rd!=1 || cd<1) wrong_arg_in("setdiag");
		for (i=0; i<r; i++)
		{	if (i+k>=0 && i+k<c)
			{	*mat(m,c,i,i+k)=*md;
				ik++; if (ik<cd) md++;
			}
		}
	}
	else if (hd->type==s_cmatrix)
	{	result=new_cmatrix(r,c,""); if (error) return;
		m=matrixof(result);
        memmove((char *)m,(char *)mhd,(long)c*r*(long)2*sizeof(double));
		getmatrix(hd2,&rd,&cd,&md);
		if (rd!=1 || cd<1) wrong_arg_in("setdiag");
		m=matrixof(result);
		for (i=0; i<r; i++)
		{	if (i+k>=0 && i+k<c) 
			{	*cmat(m,c,i,i+k)=*md;
				*(cmat(m,c,i,i+k)+1)=*(md+1);
				ik++; if (ik<cd) md+=2;
			}
		}
	}
	else wrong_arg_in("setdiag");
	moveresult(st,result);
}

void mdiag2 (header *hd)
{	header *st=hd,*hd1,*result;
	int c,r,i,n,l;
	double *m,*mh,*mr;
	hd1=next_param(hd);
	hd=getvalue(hd); if (hd1) hd1=getvalue(hd1);
	if (error) return;
	if (hd1->type!=s_real) wrong_arg();
	n=(int)*realof(hd1);
	if (hd->type==s_matrix || hd->type==s_real)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(1,r,""); if (error) return;
		mr=matrixof(result); l=0;
		for (i=0; i<r; i++)
		{	if (i+n>=c) break;
			if (i+n>=0) { l++; *mr++=*mat(m,c,i,i+n); }
		}
		dimsof(result)->c=l;
		result->size=matrixsize(1,c);
	}
	else if (hd->type==s_cmatrix || hd->type==s_complex)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(1,r,""); if (error) return;
		mr=matrixof(result); l=0;
		for (i=0; i<r; i++)
		{	if (i+n>=c) break;
			if (i+n>=0) 
			{	l++;
				mh=cmat(m,c,i,i+n);
				*mr++=*mh++;
				*mr++=*mh;
			}
		}
		dimsof(result)->c=l;
		result->size=cmatrixsize(1,c);
	}
	else wrong_arg();
	moveresult(st,result);
}

void mband (header *hd)
{	header *st=hd,*hd1,*hd2,*result;
	int i,j,c,r,n1,n2;
	double *m,*mr;
	hd1=next_param(hd);
	hd2=next_param(hd1);
	hd=getvalue(hd); hd1=getvalue(hd1); hd2=getvalue(hd2);
	if (error) return;
	if (hd1->type!=s_real || hd2->type!=s_real) wrong_arg();
	n1=(int)*realof(hd1); n2=(int)*realof(hd2);
	if (hd->type==s_matrix || hd->type==s_real)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
			for (j=0; j<c; j++)
			{	if (j-i>=n1 && j-i<=n2) *mr++=*m++;
				else { *mr++=0.0; m++; }
			}
	}
	else if (hd->type==s_cmatrix || hd->type==s_complex)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
			for (j=0; j<c; j++)
			{	if (j-i>=n1 && j-i<=n2) { *mr++=*m++; *mr++=*m++; }
				else { *mr++=0.0; *mr++=0.0; m+=2; }
			}
	}
	else wrong_arg();
	moveresult(st,result);
}

void mdup (header *hd)
{	header *result,*st=hd,*hd1;
	double x,*m,*m1,*m2;
	int c,i,n,j,r;
	hd=getvalue(hd);
	hd1=next_param(st);
	if (!hd1) wrong_arg();
	hd1=getvalue(hd1); if (error) return;
	if (hd1->type!=s_real) wrong_arg();
	x=*realof(hd1); n=(int)x;
	if (n<1 || x>=INT_MAX) wrong_arg();
	if (hd->type==s_matrix && dimsof(hd)->r==1)
	{	c=dimsof(hd)->c;
		result=new_matrix(n,c,"");
		if (error) return;
		m1=matrixof(hd); m2=matrixof(result);
		for (i=0; i<n; i++)
		{	m=mat(m2,c,i,0);
			memmove((char *)m,(char *)m1,c*sizeof(double));
		}
	}
	else if (hd->type==s_matrix && dimsof(hd)->c==1)
	{	r=dimsof(hd)->r;
		result=new_matrix(r,n,"");
		if (error) return;
		m1=matrixof(hd); m2=matrixof(result);
		for (i=0; i<r; i++)
		{	for (j=0; j<n; j++)
				*mat(m2,n,i,j)=*mat(m1,1,i,0);
		}
	}
	else if (hd->type==s_real)
	{	result=new_matrix(n,1,""); if (error) return;
		m1=matrixof(result);
		for (i=0; i<n; i++) *m1++=*realof(hd);
	}
	else if (hd->type==s_cmatrix && dimsof(hd)->r==1)
	{	c=dimsof(hd)->c;
		result=new_cmatrix(n,c,"");
		if (error) return;
		m1=matrixof(hd); m2=matrixof(result);
		for (i=0; i<n; i++)
		{	m=cmat(m2,c,i,0);
			memmove((char *)m,(char *)m1,(long)2*c*sizeof(double));
		}
	}
	else if (hd->type==s_cmatrix && dimsof(hd)->c==1)
	{	r=dimsof(hd)->r;
		result=new_cmatrix(r,n,"");
		if (error) return;
		m1=matrixof(hd); m2=matrixof(result);
		for (i=0; i<r; i++)
		{	for (j=0; j<n; j++)
				copy_complex(cmat(m2,n,i,j),cmat(m1,1,i,0));
		}
	}
	else if (hd->type==s_complex)
	{	result=new_cmatrix(n,1,""); if (error) return;
		m1=matrixof(result);
		for (i=0; i<n; i++) { *m1++=*realof(hd); *m1++=*imagof(hd); }
	}
	else if (hd->type==s_imatrix && dimsof(hd)->r==1)
	{	c=dimsof(hd)->c;
		result=new_imatrix(n,c,"");
		if (error) return;
		m1=matrixof(hd); m2=matrixof(result);
		for (i=0; i<n; i++)
		{	m=imat(m2,c,i,0);
			memmove((char *)m,(char *)m1,(long)2*c*sizeof(double));
		}
	}
	else if (hd->type==s_imatrix && dimsof(hd)->c==1)
	{	r=dimsof(hd)->r;
		result=new_imatrix(r,n,"");
		if (error) return;
		m1=matrixof(hd); m2=matrixof(result);
		for (i=0; i<r; i++)
		{	for (j=0; j<n; j++)
				copy_interval(cmat(m2,n,i,j),cmat(m1,1,i,0));
		}
	}
	else if (hd->type==s_interval)
	{	result=new_imatrix(n,1,""); if (error) return;
		m1=matrixof(result);
		for (i=0; i<n; i++) { *m1++=*aof(hd); *m1++=*bof(hd); }
	}
	else wrong_arg();
	moveresult(st,result);
}

/****************************************************************************
 *	operation on matrix rows / cols
 *
 ****************************************************************************/

void transpose (header *hd)
/***** transpose 
	transpose a matrix
*****/
{	header *hd1,*st=hd;
	double *m,*m1,*mh;
	int c,r,i,j;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		hd1=new_matrix(c,r,""); if (error) return;
		m1=matrixof(hd1);
		for (i=0; i<r; i++)
		{	mh=m1+i;
			for (j=0; j<c; j++)
			{	*mh=*m++; mh+=r;
			}
		}
	}
	else if (hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		hd1=new_cmatrix(c,r,""); if (error) return;
		m1=matrixof(hd1);
		for (i=0; i<r; i++)
		{   mh=m1+(long)2*i;
			for (j=0; j<c; j++)
			{	*mh=*m++; *(mh+1)=*m++;
                mh+=(long)2*r;
			}
		}
	}
	else if (hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		hd1=new_imatrix(c,r,""); if (error) return;
		m1=matrixof(hd1);
		for (i=0; i<r; i++)
		{   mh=m1+(long)2*i;
			for (j=0; j<c; j++)
			{	*mh=*m++; *(mh+1)=*m++;
				mh+=(long)2*r;
			}
		}
	}
	else if (hd->type==s_real || hd->type==s_complex
		|| hd->type==s_interval)
	{	hd1=hd;
	}
	else
	{	error=24; output("\' not defined for this value!\n");
		return;
	}
	moveresult(st,hd1);
}

void msum (header *hd)
{	header *st=hd,*result;
	int c,r,i,j;
	double *m,*mr,s,si,x,y;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++) 
		{	s=0.0;
			for (j=0; j<c; j++) s+=*m++;
			*mr++=s;
		}
	}
	else if (hd->type==s_complex || hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++) 
		{	s=0.0; si=0.0;
			for (j=0; j<c; j++) { s+=*m++; si+=*m++; }
			*mr++=s; *mr++=si;
		}
	}
	else if (hd->type==s_interval || hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	s=0.0; si=0.0;
			for (j=0; j<c; j++)
			{	interval_add(&s,&si,m,m+1,&x,&y);
				s=x; si=y; m+=2;
			}
			*mr++=s; *mr++=si;
		}
	}
	else wrong_arg_in("sum");
	moveresult(st,result);
}

void mprod (header *hd)
{	header *st=hd,*result;
	int c,r,i,j;
	double *m,*mr,s,si,h,hi,x,y;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++) 
		{	s=1.0;
			for (j=0; j<c; j++) s*=*m++;
			*mr++=s;
		}
	}
	else if (hd->type==s_complex || hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	s=1.0;
			for (j=0; j<c; j++)
			{	complex_multiply(&s,&si,m,m+1,&h,&hi);
				s=h; si=hi; m+=2;
			}
			*mr++=s; *mr++=si;
		}
	}
	else if (hd->type==s_interval || hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,1,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	s=1.0; si=1.0;
			for (j=0; j<c; j++)
			{	interval_mult(&s,&si,m,m+1,&x,&y);
				s=x; si=y; m+=2;
			}
			*mr++=s; *mr++=si;
		}
	}
	else wrong_arg_in("prod");
	moveresult(st,result);
}

void mcumsum (header *hd)
{	header *result,*st=hd;
	double *m,*mr,sum=0,sumr=0,sumi=0;
	int r,c,i,j;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		if (c<1) result=new_matrix(r,1,"");
		else result=new_matrix(r,c,"");
		if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++) 
		{	if (c>=1) sum=*m++;
			*mr++=sum;
			for (j=1; j<c; j++) 
			{	sum+=*m++;
				*mr++=sum;
			}
		}
	}
	else if (hd->type==s_complex || hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		if (c<1) result=new_cmatrix(r,1,"");
		else result=new_cmatrix(r,c,"");
		if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	if (c>=1) { sumr=*m++; sumi=*m++; }
			*mr++=sumr; *mr++=sumi;
			for (j=1; j<c; j++) 
			{	sumr+=*m++; *mr++=sumr;
				sumi+=*m++; *mr++=sumi;
			}
		}
	}
	else wrong_arg_in("cumsum");
	moveresult(st,result);
}

void mcumprod (header *hd)
{	header *result,*st=hd;
	double *m,*mr,sum=1,sumi=1,sumr=0,x,y;
	int r,c,i,j;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		if (c<1) result=new_matrix(r,1,"");
		else result=new_matrix(r,c,"");
		if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++) 
		{	if (c>=1) sum=*m++; 
			*mr++=sum;
			for (j=1; j<c; j++) 
			{	sum*=*m++;
				*mr++=sum;
			}
		}
	}
	else if (hd->type==s_complex || hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		if (c<1) result=new_cmatrix(r,1,"");
		else result=new_cmatrix(r,c,"");
		if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	if (c>=1) { sumr=*m++; sumi=*m++; }
			*mr++=sumr; *mr++=sumi;
			for (j=1; j<c; j++)
			{	sum=sumr*(*m)-sumi*(*(m+1));
				sumi=sumr*(*(m+1))+sumi*(*m);
				sumr=sum;
				m+=2;
				*mr++=sumr;
				*mr++=sumi;
			}
		}
	}
	else if (hd->type==s_interval || hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		if (c<1) result=new_imatrix(r,1,"");
		else result=new_imatrix(r,c,"");
		if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	if (c>=1) { sumr=*m++; sumi=*m++; }
			*mr++=sumr; *mr++=sumi;
			for (j=1; j<c; j++)
			{   interval_mult(m,m+1,&sumr,&sumi,&x,&y);
				sumr=x; sumi=y;
				m+=2;
				*mr++=sumr;
				*mr++=sumi;
			}
		}
	}
	else if (hd->type==s_interval || hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		if (c<1) result=new_imatrix(r,1,"");
		else result=new_imatrix(r,c,"");
		if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	if (c>=1) { sumr=*m++; sumi=*m++; }
			*mr++=sumr; *mr++=sumi;
			for (j=1; j<c; j++)
			{   interval_add(m,m+1,&sumr,&sumi,&x,&y);
				sumr=x; sumi=y;
				m+=2;
				*mr++=sumr;
				*mr++=sumi;
			}
		}
	}
	else wrong_arg_in("cumprod");
	moveresult(st,result);
}

void mflipx (header *hd)
{	header *st=hd,*result;
	double *m,*mr,*mr1;
	int i,j,c,r;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_complex || hd->type==s_interval)
	{	moveresult(st,hd); return;
	}
	else if (hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=mr+(c-1);
			for (j=0; j<c; j++) *mr1--=*m++;
			mr+=c;
		}
	}
	else if (hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=mr+(2l*(c-1)+1);
			for (j=0; j<c; j++)
			{	*mr1--=*(m+1); *mr1--=*m; m+=2;
			}
			mr+=2l*c;
		}
	}
	else if (hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=mr+(2l*(c-1)+1);
			for (j=0; j<c; j++)
			{	*mr1--=*(m+1); *mr1--=*m; m+=2;
			}
			mr+=2l*c;
		}
	}
	else wrong_arg_in("flipx");
	moveresult(st,result);
}

void mflipy (header *hd)
{	header *st=hd,*result;
	double *m,*mr;
	int i,c,r;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_complex || hd->type==s_interval)
	{	moveresult(st,hd); return;
	}
	else if (hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,c,""); if (error) return;
		mr=matrixof(result);
		mr+=(long)(r-1)*c;
		for (i=0; i<r; i++)
		{	memmove((char *)mr,(char *)m,c*sizeof(double));
			m+=c; mr-=c;
		}
	}
	else if (hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		mr+=2l*(long)(r-1)*c;
		for (i=0; i<r; i++)
		{	memmove((char *)mr,(char *)m,2l*c*sizeof(double));
			m+=2l*c; mr-=2l*c;
		}
	}
	else if (hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		mr+=2l*(long)(r-1)*c;
		for (i=0; i<r; i++)
		{	memmove((char *)mr,(char *)m,2l*c*sizeof(double));
			m+=2l*c; mr-=2l*c;
		}
	}
	else wrong_arg_in("flipy");
	moveresult(st,result);
}

void mrotleft (header *hd)
{	header *st=hd,*result;
	double *m,*mr,*mr1;
	int i,j,c,r;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_complex || hd->type==s_interval)
	{	moveresult(st,hd); return;
	}
	else if (hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{   mr1=m+1;
        	for (j=0; j<c-1; j++) *mr++=*mr1++;
            *mr++=*m;
            m+=c;
		}
	}
	else if (hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=m+2l;
			for (j=0; j<c-1; j++)
			{	*mr++=*mr1++; *mr++=*mr1++;
			}
            *mr++=*m; *mr++=*(m+1);
			m+=2l*c;
		}
	}
	else if (hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=m+2l;
			for (j=0; j<c-1; j++)
			{	*mr++=*mr1++; *mr++=*mr1++;
			}
            *mr++=*m; *mr++=*(m+1);
			m+=2l*c;
		}
	}
	else wrong_arg_in("flipx");
	moveresult(st,result);
}

void mrotright (header *hd)
{	header *st=hd,*result;
	double *m,*mr,*mr1;
	int i,j,c,r;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_complex || hd->type==s_interval)
	{	moveresult(st,hd); return;
	}
	else if (hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{   mr1=m;
            *mr++=*(m+c-1);
        	for (j=1; j<c; j++) *mr++=*mr1++;
            m+=c;
		}
	}
	else if (hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=m;
            *mr++=*(m+2l*(c-1)); *mr++=*(m+2l*(c-1)+1);
			for (j=1; j<c; j++)
			{	*mr++=*mr1++; *mr++=*mr1++;
			}
			m+=2l*c;
		}
	}
	else if (hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=m;
            *mr++=*(m+2l*(c-1)); *mr++=*(m+2l*(c-1)+1);
			for (j=1; j<c; j++)
			{	*mr++=*mr1++; *mr++=*mr1++;
			}
			m+=2l*c;
		}
	}
	else wrong_arg_in("flipx");
	moveresult(st,result);
}

void mshiftleft (header *hd)
{	header *st=hd,*result;
	double *m,*mr,*mr1;
	int i,j,c,r;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_complex || hd->type==s_interval)
	{	moveresult(st,hd); return;
	}
	else if (hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{   mr1=m+1;
        	for (j=0; j<c-1; j++) *mr++=*mr1++;
            *mr++=0;
            m+=c;
		}
	}
	else if (hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=m+2l;
			for (j=0; j<c-1; j++)
			{	*mr++=*mr1++; *mr++=*mr1++;
			}
            *mr++=0; *mr++=0;
			m+=2l*c;
		}
	}
	else if (hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=m+2l;
			for (j=0; j<c-1; j++)
			{	*mr++=*mr1++; *mr++=*mr1++;
			}
            *mr++=0; *mr++=0;
			m+=2l*c;
		}
	}
	else wrong_arg_in("flipx");
	moveresult(st,result);
}

void mshiftright (header *hd)
{	header *st=hd,*result;
	double *m,*mr,*mr1;
	int i,j,c,r;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_complex || hd->type==s_interval)
	{	moveresult(st,hd); return;
	}
	else if (hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{   mr1=m;
            *mr++=0;
        	for (j=1; j<c; j++) *mr++=*mr1++;
            m+=c;
		}
	}
	else if (hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_cmatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=m;
            *mr++=0; *mr++=0;
			for (j=1; j<c; j++)
			{	*mr++=*mr1++; *mr++=*mr1++;
			}
			m+=2l*c;
		}
	}
	else if (hd->type==s_imatrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_imatrix(r,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++)
		{	mr1=m;
            *mr++=0; *mr++=0;
			for (j=1; j<c; j++)
			{	*mr++=*mr1++; *mr++=*mr1++;
			}
			m+=2l*c;
		}
	}
	else wrong_arg_in("flipx");
	moveresult(st,result);
}

void mvconcat (header *hd)
{	header *st=hd,*hd1,*result;
	double *m,*m1;
	int r,c,r1,c1;
	unsigned long size=0;
	hd=getvalue(hd);
	hd1=next_param(st);
	if (!hd1) wrong_arg();
	hd1=getvalue(hd1); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	if (hd1->type==s_real || hd1->type==s_matrix)
		{	getmatrix(hd,&r,&c,&m);
			getmatrix(hd1,&r1,&c1,&m1);
			if (r!=0 && (c!=c1 || (long)r+r1>INT_MAX)) wrong_arg();
			result=new_matrix(r+r1,c1,"");
			if (r!=0)
			{	size=(long)r*c*sizeof(double);
				memmove((char *)matrixof(result),(char *)m,
					size);
			}
			memmove((char *)matrixof(result)+size,
				(char *)m1,(long)r1*c1*sizeof(double));
		}
		else if (hd1->type==s_complex || hd1->type==s_cmatrix)
		{	make_complex(st);
			mvconcat(st);
			return;
		}
		else if (hd1->type==s_interval || hd1->type==s_imatrix)
		{	make_interval(st);
			mvconcat(st);
			return;
		}
		else wrong_arg();
	}
	else if (hd->type==s_complex || hd->type==s_cmatrix)
	{	if (hd1->type==s_complex || hd1->type==s_cmatrix)
		{	getmatrix(hd,&r,&c,&m);
			getmatrix(hd1,&r1,&c1,&m1);
			if (r!=0 && (c!=c1 || (long)r+r1>INT_MAX)) wrong_arg();
			result=new_cmatrix(r+r1,c1,"");
			if (r!=0)
			{	size=(long)r*(long)2*c*sizeof(double);
				memmove((char *)matrixof(result),(char *)m,
					size);
			}
			memmove((char *)matrixof(result)+size,
				(char *)m1,(long)r1*(long)2*c1*sizeof(double));
		}
		else if (hd1->type==s_real || hd1->type==s_matrix)
		{	make_complex(next_param(st));
			mvconcat(st);
			return;
		}
		else wrong_arg();
	}
	else if (hd->type==s_interval || hd->type==s_imatrix)
	{	if (hd1->type==s_interval || hd1->type==s_imatrix)
		{	getmatrix(hd,&r,&c,&m);
			getmatrix(hd1,&r1,&c1,&m1);
			if (r!=0 && (c!=c1 || (long)r+r1>INT_MAX)) wrong_arg();
			result=new_imatrix(r+r1,c1,"");
			if (r!=0)
			{	size=(long)r*(long)2*c*sizeof(double);
				memmove((char *)matrixof(result),(char *)m,
					size);
			}
			memmove((char *)matrixof(result)+size,
				(char *)m1,(long)r1*(long)2*c1*sizeof(double));
		}
		else if (hd1->type==s_real || hd1->type==s_matrix)
		{	make_interval(next_param(st));
			mvconcat(st);
			return;
		}
		else wrong_arg();
	}
	else wrong_arg();
	moveresult(st,result);
}

void mhconcat (header *hd)
{	header *st=hd,*hd1,*result;
	double *m,*m1,*mr;
	int r,c,r1,c1,i;
	hd=getvalue(hd);
	hd1=next_param(st);
	if (!hd1) wrong_arg();
	hd1=getvalue(hd1); if (error) return;
	if (hd->type==s_string && hd1->type==s_string)
	{	result=new_string(stringof(hd),
			strlen(stringof(hd))+strlen(stringof(hd1))+1,"");
		strcat(stringof(result),stringof(hd1));
	}
	else if (hd->type==s_real || hd->type==s_matrix)
	{	if (hd1->type==s_real || hd1->type==s_matrix)
		{	getmatrix(hd,&r,&c,&m);
			getmatrix(hd1,&r1,&c1,&m1);
			if (c!=0 && (r!=r1 || (long)c+c1>INT_MAX)) wrong_arg();
			result=new_matrix(r1,c+c1,"");
			mr=matrixof(result);
			for (i=0; i<r1; i++)
			{	if (c!=0) memmove((char *)mat(mr,c+c1,i,0),
					(char *)mat(m,c,i,0),c*sizeof(double));
				memmove((char *)mat(mr,c+c1,i,c),
					(char *)mat(m1,c1,i,0),c1*sizeof(double));
			}
		}
		else if (hd1->type==s_complex || hd1->type==s_cmatrix)
		{	make_complex(st);
			mhconcat(st);
			return;
		}
		else if (hd1->type==s_interval || hd1->type==s_imatrix)
		{	make_interval(st);
			mhconcat(st);
			return;
		}
		else wrong_arg();
	}
	else if (hd->type==s_complex || hd->type==s_cmatrix)
	{	if (hd1->type==s_complex || hd1->type==s_cmatrix)
		{	getmatrix(hd,&r,&c,&m);
			getmatrix(hd1,&r1,&c1,&m1);
			if (c!=0 && (r!=r1 || (long)c+c1>INT_MAX)) wrong_arg();
			result=new_cmatrix(r1,c+c1,"");
			mr=matrixof(result);
			for (i=0; i<r1; i++)
			{	if (c!=0) memmove((char *)cmat(mr,c+c1,i,0),
					(char *)cmat(m,c,i,0),c*(long)2*sizeof(double));
				memmove((char *)cmat(mr,c+c1,i,c),
					(char *)cmat(m1,c1,i,0),c1*(long)2*sizeof(double));
			}
		}
		else if (hd1->type==s_real || hd1->type==s_matrix)
		{	make_complex(next_param(st)); if (error) return;
			mhconcat(st);
			return;
		}
		else wrong_arg();
	}
	else if (hd->type==s_interval || hd->type==s_imatrix)
	{	if (hd1->type==s_interval || hd1->type==s_imatrix)
		{	getmatrix(hd,&r,&c,&m);
			getmatrix(hd1,&r1,&c1,&m1);
			if (c!=0 && (r!=r1 || (long)c+c1>INT_MAX)) wrong_arg();
			result=new_imatrix(r1,c+c1,"");
			mr=matrixof(result);
			for (i=0; i<r1; i++)
			{	if (c!=0) memmove((char *)cmat(mr,c+c1,i,0),
					(char *)cmat(m,c,i,0),c*(long)2*sizeof(double));
				memmove((char *)cmat(mr,c+c1,i,c),
					(char *)cmat(m1,c1,i,0),c1*(long)2*sizeof(double));
			}
		}
		else if (hd1->type==s_real || hd1->type==s_matrix)
		{	make_interval(next_param(st)); if (error) return;
			mhconcat(st);
			return;
		}
		else wrong_arg();
	}
	else wrong_arg();
	moveresult(st,result);
}

/****************************************************************************
 *	comparing
 *
 ****************************************************************************/

void mscompare (header *hd)
{	header *st=hd,*hd1,*result;
	hd=getvalue(hd);
	hd1=getvalue(nextof(st));
	if (error) return;
	if (hd->type==s_string && hd1->type==s_string)
	{	result=new_real(strcmp(stringof(hd),stringof(hd1)),"");
		if (error) return;
	}
	else wrong_arg();
	moveresult(st,result);
}

void mfind (header *hd)
{	header *st=hd,*hd1,*result;
	double *m,*m1,*mr;
	int i,j,k,c,r,c1,r1;
	hd=getvalue(hd);
	hd1=getvalue(nextof(st));
	if (error) return;
	if ((hd->type!=s_matrix && hd->type!=s_real) || 
	    (hd1->type!=s_matrix && hd1->type!=s_real)) wrong_arg();
	getmatrix(hd,&c,&r,&m);
	getmatrix(hd1,&c1,&r1,&m1);
	if (c!=1 && r!=1) wrong_arg();
	if (r!=1) c=r;
	result=new_matrix(c1,r1,""); if (error) return;
	mr=matrixof(result);
	for (i=0; i<r1; i++)
		for (j=0; j<c1; j++)
		{	k=0;
			while (k<c && m[k]<=*m1) k++;
			if (k==c && *m1<=m[c-1]) k=c-1;
			*mr++=k; m1++;
		}
	moveresult(st,result);
}

void mextrema (header *hd)
{	header *result,*st=hd;
	double x,*m,*mr,min,max;
	int r,c,i,j,imin,imax;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		result=new_matrix(r,4,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<r; i++) 
		{	min=max=*m; imin=imax=0; m++;
			for (j=1; j<c; j++) 
			{	x=*m++;
				if (x<min) { min=x; imin=j; }
				if (x>max) { max=x; imax=j; }
			}
			*mr++=min; *mr++=imin+1; *mr++=max; *mr++=imax+1;
		}
	}
	else wrong_arg_in("extrema");
	moveresult(st,result);
}

void many (header *hd)
{	header *st=hd,*result;
	int c,r,res=0;
	long i,n;
	double *m;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		n=(long)(c)*r;
	}
	else if (hd->type==s_complex || hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m);
        n=(long)2*(long)(c)*r;
	}
	else wrong_arg_in("any");
	for (i=0; i<n; i++)
		if (*m++!=0.0) { res=1; break; }
	result=new_real(res,""); if (error) return;
	moveresult(st,result);
}

void mnonzeros (header *hd)
{	header *st=hd,*result;
	double *m,*mr;
	int r,c,i,k;
	hd=getvalue(hd); if (error) return;
	if (hd->type!=s_real && hd->type!=s_matrix) wrong_arg_in("nonzeros");
	getmatrix(hd,&r,&c,&m);
	if (r!=1 && c!=1) wrong_arg_in("nonzeros");
	if (c==1) c=r;
	result=new_matrix(1,c,""); if (error) return;
	k=0; mr=matrixof(result);
	for (i=0; i<c; i++)
	{	if (*m++!=0.0)
		{	*mr++=i+1; k++;
		}
	}
	dimsof(result)->c=k;
	result->size=matrixsize(1,k);
	moveresult(st,result);
}

/****************************************************************************
 *	misc
 *
 ****************************************************************************/

void cscalp (double *s, double *si, double *x, double *xi,
	double *y, double *yi);
void ccopy (double *y, double *x, double *xi);

void wmultiply (header *hd)
/***** multiply
	matrix multiplication for weakly nonzero matrices.
*****/
{	header *result,*st=hd,*hd1;
	dims *d,*d1;
	double *m,*m1,*m2,*mm1,*mm2,x,y,null=0.0;
	int i,j,c,r,k;
	hd=getvalue(hd); hd1=getvalue(nextof(st));
	if (error) return;
	if (hd->type==s_matrix && hd1->type==s_matrix)
	{	d=dimsof(hd);
		d1=dimsof(hd1);
		if (d->c != d1->r)
		{	error=8; output("Cannot multiply these!\n");
			return;
		}
		r=d->r; c=d1->c;
		result=new_matrix(r,c,"");
		if (error) return;
		m=matrixof(result);
		m1=matrixof(hd);
		m2=matrixof(hd1);
		for (i=0; i<r; i++)
			for (j=0; j<c; j++)
			{	mm1=mat(m1,d->c,i,0); mm2=m2+j;
				x=0.0;
				for (k=0; k<d->c; k++)
				{	if ((*mm1!=0.0)&&(*mm2!=0.0)) x+=(*mm1)*(*mm2);
					mm1++; mm2+=d1->c;
				}
				*mat(m,c,i,j)=x;
			}
		moveresult(st,result);
		return;
	}
	if (hd->type==s_matrix && hd1->type==s_cmatrix)
	{	d=dimsof(hd);
		d1=dimsof(hd1);
		if (d->c != d1->r)
		{	error=8; output("Cannot multiply these!\n");
			return;
		}
		r=d->r; c=d1->c;
		result=new_cmatrix(r,c,"");
		if (error) return;
		m=matrixof(result);
		m1=matrixof(hd);
		m2=matrixof(hd1);
		for (i=0; i<r; i++)
			for (j=0; j<c; j++)
            {   mm1=mat(m1,d->c,i,0); mm2=m2+(long)2*j;
				x=0.0; y=0.0;
				for (k=0; k<d->c; k++)
				{	if ((*mm2!=0.0 || *(mm2+1)!=0.0) &&
							(*mm1!=0.0))
					cscalp(&x,&y,mm1,&null,mm2,mm2+1);
					mm1++; mm2+=2*d1->c;
				}
				ccopy(cmat(m,c,i,j),&x,&y);
			}
		moveresult(st,result);
		return;
	}
	if (hd->type==s_cmatrix && hd1->type==s_matrix)
	{	d=dimsof(hd);
		d1=dimsof(hd1);
		if (d->c != d1->r)
		{	error=8; output("Cannot multiply these!\n");
			return;
		}
		r=d->r; c=d1->c;
		result=new_cmatrix(r,c,"");
		if (error) return;
		m=matrixof(result);
		m1=matrixof(hd);
		m2=matrixof(hd1);
		for (i=0; i<r; i++)
			for (j=0; j<c; j++)
			{	mm1=cmat(m1,d->c,i,0); mm2=m2+j;
				x=0.0; y=0.0;
				for (k=0; k<d->c; k++)
				{	if ((*mm2!=0.0) &&
							(*mm1!=0.0 ||  *(mm1+1)!=0.0))
					cscalp(&x,&y,mm1,mm1+1,mm2,&null);
					mm1+=2; mm2+=d1->c;
				}
				ccopy(cmat(m,c,i,j),&x,&y);
			}
		moveresult(st,result);
		return;
	}
	if (hd->type==s_cmatrix && hd1->type==s_cmatrix)
	{	d=dimsof(hd);
		d1=dimsof(hd1);
		if (d->c != d1->r)
		{	error=8; output("Cannot multiply these!\n");
			return;
		}
		r=d->r; c=d1->c;
		result=new_cmatrix(r,c,"");
		if (error) return;
		m=matrixof(result);
		m1=matrixof(hd);
		m2=matrixof(hd1);
		for (i=0; i<r; i++)
			for (j=0; j<c; j++)
            {   mm1=cmat(m1,d->c,i,0); mm2=m2+(long)2*j;
				x=0.0; y=0.0;
				for (k=0; k<d->c; k++)
				{	if ((*mm2!=0.0 || *(mm2+1)!=0.0) &&
							(*mm1!=0.0 || *(mm1+1)!=0.0))
					cscalp(&x,&y,mm1,mm1+1,mm2,mm2+1);
					mm1+=2; mm2+=2*d1->c;
				}
				ccopy(cmat(m,c,i,j),&x,&y);
			}
		moveresult(st,result);
		return;
	}
	else wrong_arg();
}

void smultiply (header *hd)
/***** multiply
	matrix multiplication for weakly nonzero symmetric matrices.
*****/
{	header *result,*st=hd,*hd1;
	dims *d,*d1;
	double *m,*m1,*m2,*mm1,*mm2,x,y,null=0.0;
	int i,j,c,r,k;
	hd=getvalue(hd); hd1=getvalue(nextof(st));
	if (error) return;
	if (hd->type==s_matrix && hd1->type==s_matrix)
	{	d=dimsof(hd);
		d1=dimsof(hd1);
		if (d->c != d1->r)
		{	error=8; output("Cannot multiply these!\n");
			return;
		}
		r=d->r; c=d1->c;
		result=new_matrix(r,c,"");
		if (error) return;
		m=matrixof(result);
		m1=matrixof(hd);
		m2=matrixof(hd1);
		for (i=0; i<r; i++)
			for (j=i; j<c; j++)
			{	mm1=mat(m1,d->c,i,0); mm2=m2+j;
				x=0.0;
				for (k=0; k<d->c; k++)
				{	if ((*mm1!=0.0)&&(*mm2!=0.0)) x+=(*mm1)*(*mm2);
					mm1++; mm2+=d1->c;
				}
				*mat(m,c,i,j)=x;
				*mat(m,c,j,i)=x;
			}
		moveresult(st,result);
		return;
	}
	if (hd->type==s_matrix && hd1->type==s_cmatrix)
	{	d=dimsof(hd);
		d1=dimsof(hd1);
		if (d->c != d1->r)
		{	error=8; output("Cannot multiply these!\n");
			return;
		}
		r=d->r; c=d1->c;
		result=new_cmatrix(r,c,"");
		if (error) return;
		m=matrixof(result);
		m1=matrixof(hd);
		m2=matrixof(hd1);
		for (i=0; i<r; i++)
			for (j=i; j<c; j++)
            {   mm1=mat(m1,d->c,i,0); mm2=m2+(long)2*j;
				x=0.0; y=0.0;
				for (k=0; k<d->c; k++)
				{	if ((*mm2!=0.0 || *(mm2+1)!=0.0) &&
							(*mm1!=0.0))
					cscalp(&x,&y,mm1,&null,mm2,mm2+1);
					mm1++; mm2+=2*d1->c;
				}
				ccopy(cmat(m,c,i,j),&x,&y); y=-y;
				ccopy(cmat(m,c,j,i),&x,&y);
			}
		moveresult(st,result);
		return;
	}
	if (hd->type==s_cmatrix && hd1->type==s_matrix)
	{	d=dimsof(hd);
		d1=dimsof(hd1);
		if (d->c != d1->r)
		{	error=8; output("Cannot multiply these!\n");
			return;
		}
		r=d->r; c=d1->c;
		result=new_cmatrix(r,c,"");
		if (error) return;
		m=matrixof(result);
		m1=matrixof(hd);
		m2=matrixof(hd1);
		for (i=0; i<r; i++)
			for (j=i; j<c; j++)
			{	mm1=cmat(m1,d->c,i,0); mm2=m2+j;
				x=0.0; y=0.0;
				for (k=0; k<d->c; k++)
				{	if ((*mm2!=0.0) &&
							(*mm1!=0.0 ||  *(mm1+1)!=0.0))
					cscalp(&x,&y,mm1,mm1+1,mm2,&null);
					mm1+=2; mm2+=d1->c;
				}
				ccopy(cmat(m,c,i,j),&x,&y); y=-y;
				ccopy(cmat(m,c,j,i),&x,&y);
			}
		moveresult(st,result);
		return;
	}
	if (hd->type==s_cmatrix && hd1->type==s_cmatrix)
	{	d=dimsof(hd);
		d1=dimsof(hd1);
		if (d->c != d1->r)
		{	error=8; output("Cannot multiply these!\n");
			return;
		}
		r=d->r; c=d1->c;
		result=new_cmatrix(r,c,"");
		if (error) return;
		m=matrixof(result);
		m1=matrixof(hd);
		m2=matrixof(hd1);
		for (i=0; i<r; i++)
			for (j=i; j<c; j++)
            {   mm1=cmat(m1,d->c,i,0); mm2=m2+(long)2*j;
				x=0.0; y=0.0;
				for (k=0; k<d->c; k++)
				{	if ((*mm2!=0.0 || *(mm2+1)!=0.0) &&
							(*mm1!=0.0 || *(mm1+1)!=0.0))
					cscalp(&x,&y,mm1,mm1+1,mm2,mm2+1);
					mm1+=2; mm2+=2*d1->c;
				}
				ccopy(cmat(m,c,i,j),&x,&y); 
				y=-y;
				ccopy(cmat(m,c,j,i),&x,&y);
			}
		moveresult(st,result);
		return;
	}
	else wrong_arg();
}