/************************************************
 Expect bugs!
 Please use and enjoy, and let me know of any bugs/mods/improvements 
 that you have found/implemented and I will fix/incorporate them into 
 this file. Thank Mr. Xushiliang once again.

					hujinshan@2002.11.3
				Airforce Engineering University
************************************************/

/*****  #include "CH4.h"  非線性方程與方程組的求解  *****/
#ifndef CH4_H_
#define CH4_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"

#include "CH3.H"
#include "CH1.H"
#include "CH13.H"
//*******************************************************************
//作參數的函數指針指向用戶自編函數，在具體應用中必須予以確定，函數定義參見文件末
//調用例子: n=ddhrt(-2.0,5.0,0.2,0.000001,x,m,ddhrtf);直接將函數名作實參即可
//*******************************************************************
int ddhrt(double a,double b,double h,double eps,double x[],int m,
		  double (*ddhrtf)(double x));//求非線性方程實根的對分法

int dnewt(double x[],double eps,int js,
		  void (*dnewtf)(double x,double y[]));//求非線性方程一個實根的牛頓法

int datkn(double x[],double eps,int js,
		  double (*datknf)(double x));//求非線性方程一個實根的埃特金迭代法

int dpqrt(double x[],double eps,
		  double (*dpqrtf)(double x));//求非線性方程一個實根的連分式解法

int dqrrt(double a[],int n,double xr[],double xi[],double eps,int jt);//求實系數代數方程全部根的QR方法
int dsrrt(double a[],int n,double xr[],double xi[]);//求實系數代數方程全部根的牛頓下山法
int dcsrt(double ar[],double ai[],int n,double xr[],double xi[]);//求復系數代數方程全部根的牛頓下山法

int dsnse(int n,double eps,double x[],int js,
		  double (*dsnsef)(double x[],double y[],int n));//求非線性組一組實根的梯度法

int dnetn(int n,double eps,double t,double h,double x[],int k,
		  void (*dnetnf)(double x[],double y[],int n));//求非線性方程一組實根的擬牛頓法

int dngin(int m,int n,double eps1,double eps2,double x[],int ka,
		  void (*dnginf)(int m,int n,double x[],double d[]),
		  void (*dngins)(int m,int n,double x[2],double* p));//求非線性方程組最小二乘解的廣義逆法

void dmtcl(double* x,double b,int m,double eps,
		   double (*dmtclf)(double x));//求非線性方程一個實根的蒙特卡洛法

void dcmtc(double* x,double* y,double b,int m,double eps,
		   double (*dcmtcf)(double x,double y));//求實函數或復函數方程一個復根的蒙特卡洛法

void dnmtc(double x[],int n,double b,int m,double eps,
		   double (*dnmtcf)(double x[],int n));//求非線性方程組一組實根的蒙特卡洛法

//*******************************************************************
int ddhrt(double a,double b,double h,double eps,double x[],int m,
		  double (*ddhrtf)(double x))//求非線性方程實根的對分法
{ //extern double ddhrtf();
  int n,js;
  double z,y,z1,y1,z0,y0;
  n=0; z=a; y=ddhrtf(z);
  while ((z<=b+h/2.0)&&(n!=m))
    { if (fabs(y)<eps)
        { n=n+1; x[n-1]=z;
          z=z+h/2.0; y=ddhrtf(z);
        }
      else
        { z1=z+h; y1=ddhrtf(z1);
          if (fabs(y1)<eps)
            { n=n+1; x[n-1]=z1;
              z=z1+h/2.0; y=ddhrtf(z);
            }
          else if (y*y1>0.0)
            { y=y1; z=z1;}
          else
            { js=0;
              while (js==0)
                { if (fabs(z1-z)<eps)
                    { n=n+1; x[n-1]=(z1+z)/2.0;
                      z=z1+h/2.0; y=ddhrtf(z);
                      js=1;
                    }
                  else
                    { z0=(z1+z)/2.0; y0=ddhrtf(z0);
                      if (fabs(y0)<eps)
                        { x[n]=z0; n=n+1; js=1;
                          z=z0+h/2.0; y=ddhrtf(z);
                        }
                      else if ((y*y0)<0.0)
                        { z1=z0; y1=y0;}
                      else { z=z0; y=y0;}
                    }
                }
            }
        }
    }
  return(n);
}
/////////////////////////////////////////////////////////////
int dnewt(double x[],double eps,int js,
		  void (*dnewtf)(double x,double y[]))//求非線性方程一個實根的牛頓法
{ //extern void dnewtf();
  int k,l;
  double y[2],d,p,x0,x1;
  l=js; k=1; x0=*x;
  dnewtf(x0,y);
  d=eps+1.0;
  while ((d>=eps)&&(l!=0))
    { if (fabs(y[1])+1.0==1.0)
        { printf("err\n"); return(-1);}
      x1=x0-y[0]/y[1];
      dnewtf(x1,y);
      d=fabs(x1-x0); p=fabs(y[0]);
      if (p>d) d=p;
      x0=x1; l=l-1;
    }
  *x=x1;
  k=js-l;
  return(k);
}
/////////////////////////////////////////////////////////////
int datkn(double x[],double eps,int js,
		  double (*datknf)(double x))//求非線性方程一個實根的埃特金迭代法
{ int flag,l;
  double u,v,x0;
  //extern double datknf();
  l=0; x0=*x; flag=0;
  while ((flag==0)&&(l!=js))
    { l=l+1; 
      u=datknf(x0); v=datknf(u);
      if (fabs(u-v)<eps) { x0=v; flag=1; }
      else x0=v-(v-u)*(v-u)/(v-2.0*u+x0);
    }
  *x=x0; l=js-l;
  return(l);
}
/////////////////////////////////////////////////////////////
int dpqrt(double x[],double eps,
		  double (*dpqrtf)(double x))//求非線性方程一個實根的連分式解法
{// extern double dpqrtf();
 int i,j,m,it,l;
 double a[10],y[10],z,h,x0,q;
 l=10; q=1.0e+35; x0=*x; h=0.0;
 while (l!=0)
   { l=l-1; j=0; it=l;
     while (j<=7)
        { if (j<=2) z=x0+0.1*j;
          else z=h;
          y[j]=dpqrtf(z);
          h=z;
          if (j==0) a[0]=z;
          else
            { m=0; i=0;
              while ((m==0)&&(i<=j-1))
                { if (fabs(h-a[i])+1.0==1.0) m=1;
                  else h=(y[j]-y[i])/(h-a[i]);
                  i=i+1;
                }
              a[j]=h;
              if (m!=0) a[j]=q;
              h=0.0;
              for (i=j-1; i>=0; i--)
                { if (fabs(a[i+1]+h)+1.0==1.0) h=q;
                  else h=-y[i]/(a[i+1]+h);
                }
              h=h+a[0];
            }
          if (fabs(y[j])>=eps) j=j+1;
          else { j=10; l=0;}
        }
      x0=h;
    }
  *x=h;
  return(10-it);
}
/////////////////////////////////////////////////////////////
int dqrrt(double a[],int n,double xr[],double xi[],double eps,int jt)
{ 
	int i,j;
	double *q;
	//extern int chhqr();
	q=(double*)malloc(n*n*sizeof(double));
	for (j=0; j<=n-1; j++)
		q[j]=-a[n-j-1]/a[n];
	for (j=n; j<=n*n-1; j++)
		q[j]=0.0;
	for (i=0; i<=n-2; i++)
		q[(i+1)*n+i]=1.0;
	i=chhqr(q,n,xr,xi,eps,jt);

	free(q); 
	return(i);
}
/////////////////////////////////////////////////////////////
static void g60(double* t,double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* p,double* q,int* k,int* it)
{ *it=1;
  while (*it==1)
    { *t=*t/1.67; *it=0;
      *x1=*x-(*t)*(*dx);
      *y1=*y-(*t)*(*dy);
      if (*k>=50)
	  { *p=sqrt((*x1)*(*x1)+(*y1)*(*y1));
          *q=exp(85.0/(*k));
          if (*p>=*q) *it=1;
        }
    }
  return;
}

static void g90(double xr[],double xi[],double a[],double *x,double* y,double* p,double* q,double* w,int* k)
//double xr[],xi[],a[];
{ int i;
  if (fabs(*y)<=1.0e-06)
    { 
	  *p=-(*x); 
	  *y=0.0; 
	  *q=0.0;
  }
  else
    { *p=-2.0*(*x); *q=(*x)*(*x)+(*y)*(*y);
      xr[*k-1]=(*x)*(*w);
      xi[*k-1]=-(*y)*(*w);
      *k=*k-1;
    }
  for (i=1; i<=*k; i++)
    { a[i]=a[i]-a[i-1]*(*p);
      a[i+1]=a[i+1]-a[i-1]*(*q);
    }
  xr[*k-1]=(*x)*(*w); xi[*k-1]=(*y)*(*w);
  *k=*k-1;
  if (*k==1)
    { xr[0]=-a[1]*(*w)/a[0]; xi[0]=0.0;}
  return;
}

static void g65(double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* dd,double* dc,double* c,int* k,int* is,int* it)
{ if (*it==0)
    { *is=1;
      *dd=sqrt((*dx)*(*dx)+(*dy)*(*dy));
      if (*dd>1.0) *dd=1.0;
      *dc=6.28/(4.5*(*k)); *c=0.0;
    }
  while(1==1)
    { *c=*c+(*dc);
      *dx=(*dd)*cos(*c); *dy=(*dd)*sin(*c);
      *x1=*x+*dx; *y1=*y+*dy;
      if (*c<=6.29)
        { *it=0; return;}
      *dd=*dd/1.67;
      if (*dd<=1.0e-07)
        { *it=1; return;}
      *c=0.0;
    }
}
//-------------------------------------------
int dsrrt(double a[],int n,double xr[],double xi[])
{ int m,i,jt,k,is,it;
  double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;
  double g,u,v,pq,g1,u1,v1;
  m=n;
  while ((m>0)&&(fabs(a[m])+1.0==1.0)) m=m-1;
  if (m<=0)
    { printf("fail\n"); return(-1);}
  for (i=0; i<=m; i++)
    a[i]=a[i]/a[m];
  for (i=0; i<=m/2; i++)
    { w=a[i]; a[i]=a[m-i]; a[m-i]=w;}
  k=m; is=0; w=1.0;
  jt=1;
  while (jt==1)
    { pq=fabs(a[k]);
	while (pq<1.0e-12)
        { xr[k-1]=0.0; xi[k-1]=0.0; k=k-1;
          if (k==1)
            { xr[0]=-a[1]*w/a[0]; xi[0]=0.0;
              return(1);
            }
          pq=fabs(a[k]);
        }
	q=log(pq); q=q/(1.0*k); q=exp(q);
      p=q; w=w*p;
      for (i=1; i<=k; i++)
        { a[i]=a[i]/q; q=q*p;}
      x=0.0001; x1=x; y=0.2; y1=y; dx=1.0;
      g=1.0e+37; 
      l40:
      u=a[0]; v=0.0;
      for (i=1; i<=k; i++)
        { p=u*x1; q=v*y1;
          pq=(u+v)*(x1+y1);
          u=p-q+a[i]; v=pq-p-q;
        }
      g1=u*u+v*v;
      if (g1>=g)
        { if (is!=0)
            { it=1;
              g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,
                  &k,&is,&it);
              if (it==0) goto l40;
            }
          else
            { g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,
                  &k,&it);
              if (t>=1.0e-03) goto l40;
              if (g>1.0e-18)
                { it=0;
                  g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,
                      &k,&is,&it);
                  if (it==0) goto l40;
                }
            }
          g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
        }
      else
        { g=g1; x=x1; y=y1; is=0;
          if (g<=1.0e-22)
	      g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
          else
            { u1=k*a[0]; v1=0.0;
              for (i=2; i<=k; i++)
                { p=u1*x; q=v1*y; pq=(u1+v1)*(x+y);
                  u1=p-q+(k-i+1)*a[i-1];
                  v1=pq-p-q;
                }
              p=u1*u1+v1*v1;
              if (p<=1.0e-20)
                { it=0;
                  g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,
                      &k,&is,&it);
                  if (it==0) goto l40;
                  g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
                }
              else
                { dx=(u*u1+v*v1)/p;
                  dy=(u1*v-v1*u)/p;
                  t=1.0+4.0/k;
                  g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,
                      &k,&it);
                  if (t>=1.0e-03) goto l40;
                  if (g>1.0e-18)
                    { it=0;
                      g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,
                          &c,&k,&is,&it);
                      if (it==0) goto l40;
                    }
                  g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
                }
            }
        }
      if (k==1) jt=0;
      else jt=1;
    }
  return(1);
}
/////////////////////////////////////////////////////////////
static void csrtg60(double* t,double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* p,double* q,int* k,int* it)
{ *it=1;
  while (*it==1)
    { *t=*t/1.67; *it=0;
      *x1=*x-*t*(*dx);
      *y1=*y-*t*(*dy);
      if (*k>=30)
	  { *p=sqrt(*x1*(*x1)+*y1*(*y1));
          *q=exp(75.0/(*k));
          if (*p>=*q) *it=1;
        }
    }
  return;
}

//-------------------------------------------
static void csrtg90(double xr[],double xi[],double ar[],double ai[],double* x,double* y,double* p,double* w,int* k)
{ 
	int i;
    for (i=1; i<=*k; i++)
      { ar[i]=ar[i]+ar[i-1]*(*x)-ai[i-1]*(*y);
        ai[i]=ai[i]+ar[i-1]*(*y)+ai[i-1]*(*x);
      }
    xr[*k-1]=*x*(*w); 
	xi[*k-1]=*y*(*w);
    *k=*k-1;
    if (*k==1)
      { *p=ar[0]*ar[0]+ai[0]*ai[0];
        xr[0]=-*w*(ar[0]*ar[1]+ai[0]*ai[1])/(*p);
        xi[0]=*w*(ar[1]*ai[0]-ar[0]*ai[1])/(*p);
      }
    return;
}
//-----------------------------------------
static void csrtg65(double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* dd,double* dc,double* c,int* k,int* is,int* it)
{ if (*it==0)
    { *is=1;
      *dd=sqrt(*dx*(*dx)+*dy*(*dy));
      if (*dd>1.0) *dd=1.0;
      *dc=6.28/(4.5*(*k)); *c=0.0;
    }
  while(1==1)
    { *c=*c+*dc;
      *dx=*dd*cos(*c); *dy=*dd*sin(*c);
      *x1=*x+*dx; *y1=*y+*dy;
      if (*c<=6.29)
        { *it=0; return;}
      *dd=*dd/1.67;
      if (*dd<=1.0e-07)
        { *it=1; return;}
      *c=0.0;
    }
}

int dcsrt(double ar[],double ai[],int n,double xr[],double xi[])
{ int m,i,jt,k,is,it;
  double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;
  double g,u,v,pq,g1,u1,v1;
  m=n;
  p=sqrt(ar[m]*ar[m]+ai[m]*ai[m]);
  while ((m>0)&&(p+1.0==1.0))
    {  m=m-1;
       p=sqrt(ar[m]*ar[m]+ai[m]*ai[m]);