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

#include "polynom.h"
#include "funcs.h"
#include "interval.h"
#include "linear.h"
#include "express.h"
#include "output.h"
#include "spread.h"
#include "mainloop.h"

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

#define max(x,y) ((x)>(y)?(x):(y))

double *polynom;
int degree;

void peval (double *x, double *r)
{	int i;
	double *p=polynom+degree,res;
	res=*p--;
	for (i=degree-1; i>=0; i--) res=res**x+(*p--);
	*r=res;
}

static void cpeval (double *x, double *xi, double *z, double *zi)
{	int i;
	double *p,h,hi;
	p=polynom+(2l*degree);
	*z=*p; *zi=*(p+1);
	p-=2;
	for (i=degree-1; i>=0; i--)
	{	complex_multiply(x,xi,z,zi,&h,&hi);
		*z= h + *p;
		*zi=hi+*(p+1); p-=2;
	}
}

void polyval (header *hd)
{	header *st=hd,*hd1,*result;
	int r,c;
	hd1=nextof(hd);
	equal_params_2(&hd,&hd1);
	getmatrix(hd,&r,&c,&polynom);
	if (r!=1) wrong_arg();
	degree=c-1;
	if (degree<0) wrong_arg();
	if (isinterval(hd)) result=map1i(ipeval,hd1);
	else result=map1(peval,cpeval,hd1);
	moveresult(st,result);
}

void polyadd (header *hd)
{	header *st=hd,*hd1,*result;
	int c,c1,c2,i,r;
	double *m1,*m2,*mr;
	complex *mc1,*mc2,*mcr;
	interval *mi1,*mi2,*mir;
	hd1=next_param(st);
	equal_params_2(&hd,&hd1); if (error) return;
	getmatrix(hd,&r,&c1,&m1); if (r!=1) wrong_arg();
	getmatrix(hd1,&r,&c2,&m2); if (r!=1) wrong_arg();
	c=max(c1,c2);
	if (iscomplex(hd)) /* complex values */
	{	mc1=(complex *)m1; mc2=(complex *)m2;
		result=new_cmatrix(1,c,""); if (error) return;
		mcr=(complex *)matrixof(result);
		for (i=0; i<c; i++)
		{	if (i>=c1) { c_copy(*mcr,*mc2); mcr++; mc2++; }
			else if (i>=c2) { c_copy(*mcr,*mc1); mcr++; mc1++; }
			else { c_add(*mc1,*mc2,*mcr); mc1++; mc2++; mcr++; }
		}
	}
	else if (isinterval(hd))
	{	mi1=(interval *)m1; mi2=(interval *)m2;
		result=new_imatrix(1,c,""); if (error) return;
		mir=(interval *)matrixof(result);
		for (i=0; i<c; i++)
		{	if (i>=c1) { i_copy(*mir,*mi2); mir++; mi2++; }
			else if (i>=c2) { i_copy(*mir,*mi1); mir++; mi1++; }
			else { i_add(*mi1,*mi2,*mir); mi1++; mi2++; mir++; }
		}
	}
	else if (isreal(hd))
	{	result=new_matrix(1,c,""); if (error) return;
		mr=matrixof(result);
		for (i=0; i<c; i++)
		{	if (i>=c1) { *mr++ = *m2++; }
			else if (i>=c2) { *mr++ = *m1++; }
			else { *mr++ = *m1++ + *m2++; }
		}
	}
	else wrong_arg();
	moveresult(st,result);
}

void polymult (header *hd)
{	header *st=hd,*hd1,*result;
	int c,c1,c2,i,r,j,k;
	double *m1,*m2,*mr,x;
	complex *mc1,*mc2,*mcr,xc,hc;
	interval *mi1,*mi2,*mir,xi,hi;
	hd1=next_param(st);
	equal_params_2(&hd,&hd1); if (error) return;
	getmatrix(hd,&r,&c1,&m1); if (r!=1) wrong_arg();
	getmatrix(hd1,&r,&c2,&m2); if (r!=1) wrong_arg();
	if ((long)c1+c2-1>INT_MAX) wrong_arg();
	c=c1+c2-1;
	if (iscomplex(hd))
	{	mc1=(complex *)m1; mc2=(complex *)m2;
		result=new_cmatrix(1,c,""); if (error) return;
		mcr=(complex *)matrixof(result);
		c_copy(xc,*mc1); mc1++;
		for (i=0; i<c2; i++) c_mult(xc,mc2[i],mcr[i]);
		for (j=1; j<c1; j++)
		{	c_copy(xc,*mc1); mc1++;
			for (k=j,i=0; i<c2-1; i++,k++)
			{	c_mult(xc,mc2[i],hc);
				c_add(hc,mcr[k],mcr[k]);
			}
			c_mult(xc,mc2[i],mcr[k]);
		}
	}
	else if (isinterval(hd))
	{	mi1=(interval *)m1; mi2=(interval *)m2;
		result=new_imatrix(1,c,""); if (error) return;
		mir=(interval *)matrixof(result);
		i_copy(xi,*mi1); mi1++;
		for (i=0; i<c2; i++) i_mult(xi,mi2[i],mir[i]);
		for (j=1; j<c1; j++)
		{	i_copy(xi,*mi1); mi1++;
			for (k=j,i=0; i<c2-1; i++,k++)
			{	i_mult(xi,mi2[i],hi);
				c_add(hi,mir[k],mir[k]);
			}
			c_mult(xi,mi2[i],mir[k]);
		}
	}
	else if (isreal(hd))
	{	result=new_matrix(1,c,""); if (error) return;
		mr=matrixof(result);
		x=*m1++;
		for (i=0; i<c2; i++) mr[i]=x*m2[i];
		for (j=1; j<c1; j++)
		{	x=*m1++;
			for (k=j,i=0; i<c2-1; i++,k++) mr[k]+=x*m2[i];
			mr[k]=x*m2[i];
		}
	}
	else wrong_arg();
	moveresult(st,result);
}

void polydiv (header *hd)
{	header *st=hd,*hd1,*result,*rest;
	int c1,c2,i,r,j;
	double *m1,*m2,*mr,*mh,x,l;
	complex *mc1,*mc2,*mcr,*mch,xc,lc,hc;
	interval *mi1,*mi2,*mir,*mih,xi,li,hi;
	hd1=next_param(st);
	equal_params_2(&hd,&hd1); if (error) return;
	getmatrix(hd,&r,&c1,&m1); if (r!=1) wrong_arg();
	getmatrix(hd1,&r,&c2,&m2); if (r!=1) wrong_arg();
	if (c1<c2)
	{	result=new_real(0.0,"");
		rest=(header *)newram;
		moveresult(rest,hd1);
	}
	else if (iscomplex(hd))
	{	mc1=(complex *)m1; mc2=(complex *)m2;
		result=new_cmatrix(1,c1-c2+1,""); if (error) return;
		mcr=(complex *)matrixof(result);
		rest=new_cmatrix(1,c2,""); if (error) return;
		mch=(complex *)newram;
		if (!freeram(c1*sizeof(complex)))
		{	output("Out of memory!\n"); error=190; return;
		}
		memmove((char *)mch,(char *)mc1,c1*sizeof(complex));
		c_copy(lc,mc2[c2-1]);
		if (lc[0]==0.0 && lc[1]==0.0) wrong_arg();
		for (i=c1-c2; i>=0; i--)
		{	c_div(mch[c2+i-1],lc,xc); c_copy(mcr[i],xc);
			for(j=0; j<c2; j++)
			{	c_mult(mc2[j],xc,hc);
				c_sub(mch[i+j],hc,mch[i+j]);
			}
		}
		memmove((char *)matrixof(rest),(char *)mch,c2*sizeof(complex));
	}
	else if (isinterval(hd))
	{	mi1=(interval *)m1; mi2=(interval *)m2;
		result=new_imatrix(1,c1-c2+1,""); if (error) return;
		mir=(interval *)matrixof(result);
		rest=new_imatrix(1,c2,""); if (error) return;
		mih=(complex *)newram;
		if (!freeram(c1*sizeof(complex)))
		{	output("Out of memory!\n"); error=190; return;
		}
		memmove((char *)mih,(char *)mi1,c1*sizeof(interval));
		i_copy(li,mi2[c2-1]);
		if (li[0]<=0.0 && li[1]>=0.0) wrong_arg();
		for (i=c1-c2; i>=0; i--)
		{	i_div(mih[c2+i-1],li,xi); c_copy(mir[i],xi);
			for(j=0; j<c2; j++)
			{	i_mult(mi2[j],xi,hi);
				i_sub(mih[i+j],hi,mih[i+j]);
			}
		}
		memmove((char *)matrixof(rest),(char *)mih,c2*sizeof(interval));
	}
	else if (isreal(hd))
	{	result=new_matrix(1,c1-c2+1,""); if (error) return;
		mr=matrixof(result);
		rest=new_matrix(1,c2,""); if (error) return;
		mh=(double *)newram;
		if (!freeram(c1*sizeof(double)))
		{	output("Out of memory!\n"); error=190; return;
		}
		memmove((char *)mh,(char *)m1,c1*sizeof(double));
		l=m2[c2-1];
		if (l==0.0) wrong_arg();
		for (i=c1-c2; i>=0; i--)
		{	x=mh[c2+i-1]/l; mr[i]=x;
			for(j=0; j<c2; j++) mh[i+j]-=m2[j]*x;
		}
		memmove((char *)matrixof(rest),(char *)mh,c2*sizeof(double));
	}
	else wrong_arg();
	moveresult(st,result);
	moveresult(nextof(st),rest);
}

void dd (header *hd)
{	header *st=hd,*hd1,*result;
	int c1,c2,i,j,r;
	double *m1,*m2,*mr;
	complex *mc1,*mc2,*mcr,hc1,hc2;
	interval *mi1,*mi2,*mir,hi1,hi2;
	hd1=next_param(st);
	equal_params_2(&hd,&hd1); if (error) return;
	getmatrix(hd,&r,&c1,&m1); if (r!=1) wrong_arg();
	getmatrix(hd1,&r,&c2,&m2); if (r!=1) wrong_arg();
	if (c1!=c2) wrong_arg();
	if (iscomplex(hd)) /* complex values */
	{	mc1=(complex *)m1; mc2=(complex *)m2;
		result=new_cmatrix(1,c1,""); if (error) return;
		mcr=(complex *)matrixof(result);
		memmove((char *)mcr,(char *)mc2,c1*sizeof(complex));
		for (i=1; i<c1; i++)
		{	for (j=c1-1; j>=i; j--)
			{	if (mc1[j][0]==mc1[j-i][0] &&
					mc1[j][1]==mc1[j-i][1]) wrong_arg();
				c_sub(mcr[j],mcr[j-1],hc1);
				c_sub(mc1[j],mc1[j-i],hc2);
				c_div(hc1,hc2,mcr[j]);
			}
		}
	}
	else if (isinterval(hd)) /* complex values */
	{	mi1=(complex *)m1; mi2=(complex *)m2;
		result=new_imatrix(1,c1,""); if (error) return;
		mir=(interval *)matrixof(result);
		memmove((char *)mir,(char *)mi2,c1*sizeof(interval));
		for (i=1; i<c1; i++)
		{	for (j=c1-1; j>=i; j--)
			{	i_sub(mir[j],mir[j-1],hi1);
				if (hi1[0]<=0 && hi1[1]>=0)
				{	output("Interval points coincide\n");
					error=1; return;
				}
				i_sub(mi1[j],mi1[j-i],hi2);
				i_div(hi1,hi2,mir[j]);
			}
		}
	}
	else if (isreal(hd))
	{	result=new_matrix(1,c1,""); if (error) return;
		mr=matrixof(result);
		memmove((char *)mr,(char *)m2,c1*sizeof(double));
		for (i=1; i<c1; i++)
		{	for (j=c1-1; j>=i; j--)
			{	if (m1[j]==m1[j-i]) wrong_arg();
				mr[j]=(mr[j]-mr[j-1])/(m1[j]-m1[j-i]);
			}
		}	
	}
	else wrong_arg();
	moveresult(st,result);
}

double *divx,*divdif;

static void rddeval (double *x, double *r)
{	int i;
	double *p=divdif+degree,res;
	res=*p--;
	for (i=degree-1; i>=0; i--) res=res*(*x-divx[i])+(*p--);
	*r=res;
}

static void cddeval (double *x, double *xi, double *z, double *zi)
{	int i;
	double *p,h,hi,*dd,xh,xhi;
	p=divdif+(2l*degree);
	dd=divx+(2l*(degree-1));
	*z=*p; *zi=*(p+1);
	p-=2;
	for (i=degree-1; i>=0; i--)
	{	xh=*x-*dd;
		xhi=*xi-*(dd+1); dd-=2;
		complex_multiply(&xh,&xhi,z,zi,&h,&hi);
		*z= h + *p;
		*zi=hi+*(p+1); p-=2;
	}
}

void ddval (header *hd)
{	header *st=hd,*hdd,*hd1,*result;
	int r,c,cd;
	hdd=nextof(st);
	hd1=nextof(hdd);
	equal_params_3(&hd,&hdd,&hd1); if (error) return;
	getmatrix(hd,&r,&c,&divx); if (r!=1 || c<1) wrong_arg();
	getmatrix(hdd,&r,&cd,&divdif); if (r!=1 || c!=cd) wrong_arg();
	degree=c-1;
	if (isinterval(hd)) result=map1i(iddeval,hd1);
	else result=map1(rddeval,cddeval,hd1);
	if (error) return;
	moveresult(st,result);
}

void polyzeros (header *hd)
{	header *st=hd,*result;
	int i,j,r,c;
	double *m,*mr,x;
	complex *mc,*mcr,xc,hc;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_real || hd->type==s_matrix)
	{	getmatrix(hd,&r,&c,&m);
		if (r!=1) wrong_arg();
		result=new_matrix(1,c+1,""); if (error) return;
		mr=matrixof(result);
		mr[0]=-m[0]; mr[1]=1.0;
		for (i=1; i<c; i++)
		{	x=-m[i]; mr[i+1]=1.0;
			for (j=i; j>=1; j--) mr[j]=mr[j-1]+x*mr[j];
			mr[0]*=x;
		}
	}
	else if (hd->type==s_complex || hd->type==s_cmatrix)
	{	getmatrix(hd,&r,&c,&m); mc=(complex *)m;
		if (r!=1) wrong_arg();
		result=new_cmatrix(1,c+1,""); if (error) return;
		mcr=(complex *)matrixof(result);
		mcr[0][0]=-mc[0][0]; mcr[0][1]=-mc[0][1];
		mcr[1][0]=1.0; mcr[1][1]=0.0;
		for (i=1; i<c; i++)
		{	xc[0]=-mc[i][0]; xc[1]=-mc[i][1];
			mcr[i+1][0]=1.0; mcr[i+1][1]=0.0;
			for (j=i; j>=1; j--) 
			{	c_mult(xc,mcr[j],hc);
				c_add(hc,mcr[j-1],mcr[j]);
			}
			c_mult(xc,mcr[0],mcr[0]);
		}
	}
	else wrong_arg();
	moveresult(st,result);
}

void polydd (header *hd)
{	header *st=hd,*hd1,*result;
	int c1,c2,i,j,r;
	double *m1,*m2,*mr,x;
	complex *mc1,*mc2,*mcr,hc,xc;
	hd1=next_param(st);
	equal_params_2(&hd,&hd1); if (error) return;
	getmatrix(hd,&r,&c1,&m1); if (r!=1) wrong_arg();
	getmatrix(hd1,&r,&c2,&m2); if (r!=1) wrong_arg();
	if (c1!=c2) wrong_arg();
	if (iscomplex(hd)) /* complex values */
	{	mc1=(complex *)m1; mc2=(complex *)m2;
		result=new_cmatrix(1,c1,""); if (error) return;
		mcr=(complex *)matrixof(result);
		c_copy(mcr[c1-1],mc2[c1-1]);
		for (i=c1-2; i>=0; i--)
		{	c_copy(xc,mc1[i]);
			c_mult(xc,mcr[i+1],hc);
			c_sub(mc2[i],hc,mcr[i]);
			for (j=i+1; j<c1-1; j++) 
			{	c_mult(xc,mcr[j+1],hc);
				c_sub(mcr[j],hc,mcr[j]);
			}
		}
	}
	else
	{	result=new_matrix(1,c1,""); if (error) return;
		mr=matrixof(result);
		mr[c1-1]=m2[c1-1];
		for (i=c1-2; i>=0; i--)
		{	x=m1[i];
			mr[i]=m2[i]-x*mr[i+1];
			for (j=i+1; j<c1-1; j++) mr[j]=mr[j]-x*mr[j+1];
		}
	}
	moveresult(st,result);
}

void polytrunc (header *hd)
{	header *st=hd,*result;
	double *m;
	complex *mc;
	int i;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_matrix && dimsof(hd)->r==1)
	{	m=matrixof(hd);
		for (i=dimsof(hd)->c-1; i>=0; i--)
		{	if (fabs(m[i])>epsilon) break;
		}
		if (i<0) result=new_real(0.0,"");
		else 
		{	result=new_matrix(1,i+1,"");
			memmove((char *)matrixof(result),(char *)matrixof(hd),
				(i+1)*sizeof(double));
		}
	}
	else if (hd->type==s_complex && dimsof(hd)->r==1)
	{	mc=(complex *)matrixof(hd);
		for (i=dimsof(hd)->c-1; i>=0; i--)
		{	if (fabs(mc[i][0])>epsilon && fabs(mc[i][1])>epsilon) 
				break;
		}
		if (i<0) result=new_complex(0.0,0.0,"");
		else 
		{	result=new_cmatrix(1,i+1,"");
			memmove((char *)matrixof(result),(char *)matrixof(hd),
				(i+1)*sizeof(complex));
		}
	}
	else wrong_arg();
	moveresult(st,result);
}

/************** bauhuber algorithm ***************/

#define ITERMAX 200
#define EPS (64*DBL_EPSILON)
#define QR 0.1
#define QI 0.8
#define EPSROOT (64*epsilon)
#define BETA (2096*EPSROOT)

static char *ram;

static void quadloes (double ar, double ai, double br, double bi,
	double cr, double ci, double *treal, double *timag)
{	double pr,pi,qr,qi,h;
	pr=br*br-bi*bi; pi=2*br*bi;
	qr=ar*cr-ai*ci; qi=ar*ci+ai*cr;
	pr=pr-4*qr; pi=pi-4*qi;
	h=sqrt(pr*pr+pi*pi);
	qr=h+pr; if (qr<0.0) qr=0; 
	qr=sqrt(qr/2);
	qi=h-pr; if (qi<0.0) qi=0; 
	qi=sqrt(qi/2);
	if (pi<0.0) qi=-qi;
	h=qr*br+qi*bi;
	if (h>0.0) { qr=-qr; qi=-qi; }
	pr=qr-br; pi=qi-bi;
	h=pr*pr+pi*pi;
	*treal=2*(cr*pr+ci*pi)/h;
	*timag=2*(ci*pr-cr*pi)/h;
}

static int cxdiv (double ar, double ai, double br, double bi,
	double *cr, double *ci)
{	double temp;
	if (br==0.0 && bi==0.0) return 1;
	if (fabs(br)>fabs(bi))
	{	temp=bi/br; br=temp*bi+br;
		*cr=(ar+temp*ai)/br;
		*ci=(ai-temp*ar)/br;
	}
	else
	{	temp=br/bi; bi=temp*br+bi;
		*cr=(temp*ar+ai)/bi;
		*ci=(temp*ai-ar)/bi;
	}
	return 0;
}

static double cxxabs (double ar, double ai)
{	if (ar==0.0) return fabs(ai);
	if (ai==0.0) return fabs(ar);
	return sqrt(ai*ai+ar*ar);
}

static void chorner (int n, int iu, double *ar, double *ai,
	double xr, double xi, double *pr, double *pi,
	double *p1r, double *p1i, double *p2r, double *p2i,
	double *rf1)
{	register int i,j;
	int i1;
	double temp,hh,tempr=0.0,tempi=0.0;
	*pr=ar[n]; *pi=ai[n];
	*p1r=*p2r=0.0; *p1i=*p2i=0.0;
	*rf1=cxxabs(*pr,*pi);
	i1=n-iu;
	for (j=n-iu,i=n-1; i>=iu; i--,j--)
	{	if (i<n-1)
		{	tempr=*p1r; tempi=*p1i;
			*p1r=*p1r * xr - *p1i * xi;
			*p1i=*p1i * xr + tempr * xi;
		}
		*p1r+=*pr; *p1i+=*pi;
		temp=*pr;
		*pr=*pr * xr - *pi * xi + ar[i];
		*pi=*pi * xr + temp * xi + ai[i];
		temp=cxxabs(*p1r,*p1i);
		hh=cxxabs(*pr,*pi); if (hh>temp) temp=hh;
		if (temp>*rf1)
		{	*rf1=temp; i1=j-1;
		}
		if (i<n-1)
		{	temp=*p2r;
			*p2r=*p2r * xr - *p2i * xi + tempr*2;
			*p2i=*p2i * xr + temp * xi + tempi*2;
		}
	}
	temp=cxxabs(xr,xi);
	if (temp!=0.0)
		*rf1=pow(temp,(double)i1)*(i1+1);
	else
		*rf1=cxxabs(*p1r,*p1i);
	*rf1*=EPS;
}

#if 0
static void scpoly (int n, double *ar, double *ai, double *scal)
{	double p,h;
	int i;
	*scal=0.0;
	p=cxxabs(ar[n],ai[n]);
	for (i=0; i<n; i++)
	{	ai[i]/=p; ar[i]/=p;
		h=pow(cxxabs(ar[i],ai[i]),1.0/(n-i));
		if (h>*scal) *scal=h;
	}
	ar[n]/=p; ai[n]/=p;
	if (*scal==0.0) *scal=1.0;
	for (p=1.0,i=n-1; i>=0; i--)
	{	p*= *scal;
		ar[i]/=p; ai[i]/=p;
	}
}
#endif

static void bauroot (int n, int iu, double *ar, double *ai, double *x0r,
	double *x0i)
{	int iter=0,i=0,aborted=0;
	double xoldr,xoldi,xnewr,xnewi,h,h1,h2,h3,h4,dzmax,dzmin,
		dxr=1,dxi=0,tempr,tempi,abs_pold,abs_pnew,abs_p1new,
		temp,ss,u,v,
		pr,pi,p1r,p1i,p2r,p2i,abs_pnoted=-1;
		
	dxr=dxi=xoldr=xoldi=0.0;
	if (n-iu==1)
	{	quadloes(0.0,0.0,ar[n],ai[n],
			ar[n-1],ai[n-1],x0r,x0i);
		goto stop;
	}
	if (n-iu==2)
	{	quadloes(ar[n],ai[n],ar[n-1],ai[n-1],
			ar[n-2],ai[n-2],x0r,x0i);
		goto stop;
	}
	xnewr=*x0r; xnewi=*x0i;
	chorner(n,iu,ar,ai,xnewr,xnewi,&pr,&pi,&p1r,&p1i,&p2r,&p2i,&ss);
	iter++;
	abs_pnew=cxxabs(pr,pi);
	if (abs_pnew==0) goto stop;
	abs_pold=abs_pnew;
	dzmin=BETA*(1+cxxabs(xnewr,xnewi));
	while (!aborted)
	{	abs_p1new=cxxabs(p1r,p1i);
		iter++;
		if (abs_pnew>abs_pold) /* Spiraling */
		{	i=0;
			temp=dxr;
			dxr=QR*dxr-QI*dxi;
			dxi=QR*dxi+QI*temp;
		}
		else /* Newton step */
		{	
			dzmax=1.0+cxxabs(xnewr,xnewi);
			h1=p1r*p1r-p1i*p1i-pr*p2r+pi*p2i;
			h2=2*p1r*p1i-pr*p2i-pi*p2r;
			if (abs_p1new>10*ss && cxxabs(h1,h2)>100*ss*ss)
				/* do a Newton step */
			{	i++;
				if (i>2) i=2;
				tempr=pr*p1r-pi*p1i;
				tempi=pr*p1i+pi*p1r;
				cxdiv(-tempr,-tempi,h1,h2,&dxr,&dxi);
				if (cxxabs(dxr,dxi)>dzmax)
				{	temp=dzmax/cxxabs(dxr,dxi);
					dxr*=temp; dxi*=temp;
					i=0;
				}
				if (i==2 && cxxabs(dxr,dxi)<dzmin/EPSROOT &&
					cxxabs(dxr,dxi)>0)
				{	i=0;
					cxdiv(xnewr-xoldr,xnewi-xoldi,dxr,dxi,&h3,&h4);
					h3+=1;
					h1=h3*h3-h4*h4;
					h2=2*h3*h4;
					cxdiv(dxr,dxi,h1,h2,&h3,&h4);
					if (cxxabs(h3,h4)<50*dzmin)
					{	dxr+=h3; dxi+=h4;
					}
				}
				xoldr=xnewr; xoldi=xnewi;
				abs_pold=abs_pnew;
			}
			else /* saddle point, minimize into direction pr+i*pi */
			{	i=0;
				h=dzmax/abs_pnew;
				dxr=h*pr; dxi=h*pi;
				xoldr=xnewr; xoldi=xnewi;
				abs_pold=abs_pnew;
				do
				{	chorner(n,iu,ar,ai,xnewr+dxr,xnewi+dxi,&u,&v,
						&h,&h1,&h2,&h3,&h4);
					dxr*=2; dxi*=2;
				}
				while (fabs(cxxabs(u,v)/abs_pnew-1)<EPSROOT);
			}
		} /* end of Newton step */
		xnewr=xoldr+dxr;
		xnewi=xoldi+dxi;
		dzmin=BETA*(1+cxxabs(xoldr,xoldi));
		chorner(n,iu,ar,ai,xnewr,xnewi,&pr,&pi,
			&p1r,&p1i,&p2r,&p2i,&ss);
		abs_pnew=cxxabs(pr,pi);
		if (abs_pnew==0.0) break;
		if (cxxabs(dxr,dxi)<dzmin && abs_pnew<1e-5
			&& iter>5) break;
		if (iter>ITERMAX)
		{	iter=0;
			if (abs_pnew<=abs_pnoted) break;
			abs_pnoted=abs_pnew;
			if (test_key()==escape) { error=700; return; }
		}
	}
	*x0r=xnewr; *x0i=xnewi;
	stop: ;
/*
	chorner(n,iu,ar,ai,*x0r,*x0i,&pr,&pi,&p1r,&p1i,&p2r,&p2i,&ss);
	abs_pnew=cxxabs(pr,pi);
	printf("%20.5e +i* %20.5e, %20.5e\n",
		*x0r,*x0i,abs_pnew);
*/
}

static void _polydiv (int n, int iu, double *ar, double *ai,
	double x0r, double x0i)
{	int i;
	for (i=n-1; i>iu; i--)
	{	ar[i]+=ar[i+1]*x0r-ai[i+1]*x0i;
		ai[i]+=ai[i+1]*x0r+ar[i+1]*x0i;
	}
}

static void bauhuber (double *p, int n, double *result, int all,
	double startr, double starti)
{	double *ar,*ai,scalefak=1.0;
	int i;
	double x0r,x0i;
	ram=newram;
	if (!freeram(2*(n+1)*sizeof(double))) outofram();
	ar=(double *)ram;
	ai=ar+n+1;
	for (i=0; i<=n; i++)
	{	ar[i]=p[2*i];
		ai[i]=p[2*i+1];
	}
/*	scpoly(n,ar,ai,&scalefak); */
	/* scalefak=1; */
	x0r=startr; x0i=starti;
	for (i=0; i<(all?n:1); i++)
	{	bauroot(n,i,ar,ai,&x0r,&x0i);
		ar[i]=scalefak*x0r;
		ai[i]=scalefak*x0i;
		if (error) 
		{	output("Bauhuber-Iteration failed!\n"); 
			error=311; return;
		}
		_polydiv(n,i,ar,ai,x0r,x0i);
		x0i=-x0i;
	}
	for (i=0; i<n; i++)
	{	result[2*i]=ar[i]; result[2*i+1]=ai[i];
	}
}

void mzeros (header *hd)
{	header *st=hd,*result;
	int r,c;
	double *m;
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_matrix) 
	{	make_complex(st); if (error) return;
		hd=getvalue(st); if (error) return; 
	}
	if (hd->type!=s_cmatrix || dimsof(hd)->r!=1 || dimsof(hd)->c<2)
	{	output("Need a complex polynomial\n"); error=300; return; }
	getmatrix(hd,&r,&c,&m);
	result=new_cmatrix(1,c-1,""); if (error) return;
	bauhuber(m,c-1,matrixof(result),1,0,0);
	moveresult(st,result);
}

void mzeros1 (header *hd)
{	header *st=hd,*hd1,*result;
	int r,c;
	double *m,xr,xi;
	hd1=nextof(hd);
	hd=getvalue(hd); if (error) return;
	if (hd->type==s_matrix)
	{	make_complex(st); if (error) return;
		hd=getvalue(st); if (error) return;
	}
	hd1=getvalue(hd1); if (error) return;
	if (hd1->type==s_real)
	{	xr=*realof(hd1); xi=0;
	}
	else if (hd1->type==s_complex)
	{	xr=*realof(hd1); xi=*(realof(hd1)+1);
	}
	else
	{	output("Need a starting value!\n"); error=300; return; }
	if (hd->type!=s_cmatrix || dimsof(hd)->r!=1 || dimsof(hd)->c<2)
	{	output("Need a complex polynomial\n"); error=300; return; }
	getmatrix(hd,&r,&c,&m);
	result=new_complex(0,0,""); if (error) return;
	bauhuber(m,c-1,realof(result),0,xr,xi);
	moveresult(st,result);
}