}
  if (m<=0)
    { printf("fail\n"); return(-1);}
  for (i=0; i<=m; i++)
    { ar[i]=ar[i]/p; ai[i]=ai[i]/p;}
  for (i=0; i<=m/2; i++)
    { w=ar[i]; ar[i]=ar[m-i]; ar[m-i]=w;
      w=ai[i]; ai[i]=ai[m-i]; ai[m-i]=w;
    }
  k=m; is=0; w=1.0;
  jt=1;
  while (jt==1)
    { pq=sqrt(ar[k]*ar[k]+ai[k]*ai[k]);
	while (pq<1.0e-12)
        { xr[k-1]=0.0; xi[k-1]=0.0; 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(1);
            }
          pq=sqrt(ar[k]*ar[k]+ai[k]*ai[k]);
        }
	q=log(pq); q=q/(1.0*k); q=exp(q);
      p=q; w=w*p;
      for (i=1; i<=k; i++)
        { ar[i]=ar[i]/q; ai[i]=ai[i]/q; q=q*p;}
      x=0.0001; x1=x; y=0.2; y1=y; dx=1.0;
      g=1.0e+37; 
      l40:
      u=ar[0]; v=ai[0];
      for (i=1; i<=k; i++)
        { p=u*x1; q=v*y1;
          pq=(u+v)*(x1+y1);
          u=p-q+ar[i]; v=pq-p-q+ai[i];
        }
      g1=u*u+v*v;
      if (g1>=g)
        { if (is!=0)
            { it=1;
              csrtg65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,
		    &c,&k,&is,&it);
              if (it==0) goto l40;
            }
          else
            { csrtg60(&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;
                  csrtg65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,
			&c,&k,&is,&it);
                  if (it==0) goto l40;
                }
            }
          csrtg90(xr,xi,ar,ai,&x,&y,&p,&w,&k);
        }
      else
        { g=g1; x=x1; y=y1; is=0;
          if (g<=1.0e-22)
	      csrtg90(xr,xi,ar,ai,&x,&y,&p,&w,&k);
          else
            { u1=k*ar[0]; v1=ai[0];
              for (i=2; i<=k; i++)
                { p=u1*x; q=v1*y; pq=(u1+v1)*(x+y);
                  u1=p-q+(k-i+1)*ar[i-1];
                  v1=pq-p-q+(k-i+1)*ai[i-1];
                }
              p=u1*u1+v1*v1;
              if (p<=1.0e-20)
                { it=0;
                  csrtg65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,
                      &c,&k,&is,&it);
                  if (it==0) goto l40;
                  csrtg90(xr,xi,ar,ai,&x,&y,&p,&w,&k);
                }
              else
                { dx=(u*u1+v*v1)/p;
                  dy=(u1*v-v1*u)/p;
                  t=1.0+4.0/k;
                  csrtg60(&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;
                      csrtg65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,
                          &c,&k,&is,&it);
                      if (it==0) goto l40;
                    }
                  csrtg90(xr,xi,ar,ai,&x,&y,&p,&w,&k);
                }
            }
        }
      if (k==1) jt=0;
      else jt=1;
    }
  return(1);
}
/////////////////////////////////////////////////////////////
int dsnse(int n,double eps,double x[],int js,
		  double (*dsnsef)(double x[],double y[],int n))//求非線性組一組實根的梯度法
{ //extern double dsnsef();
  int l,j;
  double f,d,s,*y;
  y=(double *)malloc(n*sizeof(double));
  l=js;
  f=dsnsef(x,y,n);
  while (f>=eps)
    { l=l-1;
      if (l==0) { free(y); return(js);}
      d=0.0;
      for (j=0; j<=n-1; j++) d=d+y[j]*y[j];
      if (d+1.0==1.0) { free(y); return(-1);}
      s=f/d;
      for (j=0; j<=n-1; j++) x[j]=x[j]-s*y[j];
      f=dsnsef(x,y,n);
    }
  free(y); return(js-l);
}
/////////////////////////////////////////////////////////////
int dnetn(int n,double eps,double t,double h,double x[],int k,
		  void (*dnetnf)(double x[],double y[],int n))//求非線性方程一組實根的擬牛頓法
{ 
	//extern void dnetnf();
    //extern int agaus();
    int i,j,l;
    double am,z,beta,d,*y,*a,*b;
    y=(double*)malloc(n*sizeof(double));
    a=(double*)malloc(n*n*sizeof(double));
    b=(double*)malloc(n*sizeof(double));
    l=k; am=1.0+eps;
    while (am>=eps)
      { dnetnf(x,b,n);
        am=0.0;
        for (i=0; i<=n-1; i++)
          { z=fabs(b[i]);
            if (z>am) am=z;
          }
        if (am>=eps)
          { l=l-1;
            if (l==0)
              { free(y); free(b); free(a);
                printf("fail\n"); return(0);
              }
            for (j=0; j<=n-1; j++)
              { z=x[j]; x[j]=x[j]+h;
                dnetnf(x,y,n);
                for (i=0; i<=n-1; i++) a[i*n+j]=y[i];
                x[j]=z;
              }
            if (agaus(a,b,n)==0)
              { free(y); free(a); free(b); return(-1);}
            beta=1.0;
            for (i=0; i<=n-1; i++) beta=beta-b[i];
            if (fabs(beta)+1.0==1.0)
              { free(y); free(a); free(b);
                printf("fail\n"); return(-2);
              }
            d=h/beta;
            for (i=0; i<=n-1; i++) x[i]=x[i]-d*b[i];
            h=t*h;
          }
      }
    free(y); free(a); free(b);
    return(k-l);
}
/////////////////////////////////////////////////////////////
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))//求非線性方程組最小二乘解的廣義逆法
{ 
	//extern void dnginf();
    //extern void dngins();
    //extern int agmiv();
    int i,j,k,l,kk,jt;
    double y[10],b[10],alpha,z,h2,y1,y2,y3,y0,h1;
    double *p,*d,*pp,*dx,*u,*v,*w;
    p=(double*)malloc(m*n*sizeof(double));
    d=(double*)malloc(m*sizeof(double));
    pp=(double*)malloc(n*m*sizeof(double));
    dx=(double*)malloc(n*sizeof(double));
    u=(double*)malloc(m*m*sizeof(double));
    v=(double*)malloc(n*n*sizeof(double));
    w=(double*)malloc(ka*sizeof(double));
    l=60; alpha=1.0;
    while (l>0)
      { dnginf(m,n,x,d);
        dngins(m,n,x,p);
        jt=agmiv(p,m,n,d,dx,pp,eps2,u,v,ka);
        if (jt<0)
          { free(p); free(d); free(pp); free(w);
            free(dx); free(u); free(v); return(jt);
          }
        j=0; jt=1; h2=0.0;
        while (jt==1)
          { jt=0;
            if (j<=2) z=alpha+0.01*j;
            else z=h2;
            for (i=0; i<=n-1; i++) w[i]=x[i]-z*dx[i];
            dnginf(m,n,w,d);
            y1=0.0;
            for (i=0; i<=m-1; i++) y1=y1+d[i]*d[i];
            for (i=0; i<=n-1; i++)
              w[i]=x[i]-(z+0.00001)*dx[i];
            dnginf(m,n,w,d);
            y2=0.0;
            for (i=0; i<=m-1; i++) y2=y2+d[i]*d[i];
            y0=(y2-y1)/0.00001;
            if (fabs(y0)>1.0e-10)
              { h1=y0; h2=z;
                if (j==0) { y[0]=h1; b[0]=h2;}
                else
                  { y[j]=h1; kk=0; k=0;
                    while ((kk==0)&&(k<=j-1))
                        { y3=h2-b[k];
                          if (fabs(y3)+1.0==1.0) kk=1;
                          else h2=(h1-y[k])/y3;
                          k=k+1;
                        }
                    b[j]=h2;
                    if (kk!=0) b[j]=1.0e+35;
                    h2=0.0;
                    for (k=j-1; k>=0; k--)
                      h2=-y[k]/(b[k+1]+h2);
                    h2=h2+b[0];
                  }
                j=j+1;
                if (j<=7) jt=1;
                else z=h2;
              }
          }
        alpha=z; y1=0.0; y2=0.0;
        for (i=0; i<=n-1; i++)
          { dx[i]=-alpha*dx[i];
            x[i]=x[i]+dx[i];
            y1=y1+fabs(dx[i]);
            y2=y2+fabs(x[i]);
          }
        if (y1<eps1*y2)
          { free(p); free(pp); free(d); free(w);
            free(dx); free(u); free(v); return(1);
          }
        l=l-1;
      }
    free(p); free(pp); free(d); free(dx);
    free(u); free(v); free(w); return(0);
}
/////////////////////////////////////////////////////////////
void dmtcl(double* x,double b,int m,double eps,
		   double (*dmtclf)(double x))//求非線性方程一個實根的蒙特卡洛法
{ //extern double mrnd1();
  //  extern double dmtclf();
    int k;
    double xx,a,r,y,x1,y1;
    a=b; k=1; r=1.0; xx=*x; y=dmtclf(xx);
    while (a>=eps)
      { x1=mrnd1(&r); x1=-a+2.0*a*x1;
        x1=xx+x1; y1=dmtclf(x1);
        k=k+1;
        if (fabs(y1)>=fabs(y))
          { if (k>m) { k=1; a=a/2.0; }}
        else
          { k=1; xx=x1; y=y1;
            if (fabs(y)<eps)
              { *x=xx; return; }
          }
      }
    *x=xx; return;
}
/////////////////////////////////////////////////////////////
void dcmtc(double* x,double* y,double b,int m,double eps,
		   double (*dcmtcf)(double x,double y))//求實函數或復函數方程一個復根的蒙特卡洛法
{ //extern double mrnd1();
    //extern double dcmtcf();
    int k;
    double xx,yy,a,r,z,x1,y1,z1;
    a=b; k=1; r=1.0; xx=*x; yy=*y;
    z=dcmtcf(xx,yy);
    while (a>=eps)
      { x1=-a+2.0*a*mrnd1(&r); x1=xx+x1; 
        y1=-a+2.0*a*mrnd1(&r); y1=yy+y1;
        z1=dcmtcf(x1,y1);
        k=k+1;
        if (z1>=z)
          { if (k>m) { k=1; a=a/2.0; }}
        else
          { k=1; xx=x1; yy=y1;  z=z1;
            if (z<eps)
              { *x=xx; *y=yy; return; }
          }
      }
    *x=xx; *y=yy; return;
}
/////////////////////////////////////////////////////////////
void dnmtc(double x[],int n,double b,int m,double eps,
		   double (*dnmtcf)(double x[],int n))//求非線性方程組一組實根的蒙特卡洛法
{ 
	//extern double mrnd1();
    //extern double dnmtcf();
    int k,i;
    double a,r,*y,z,z1;
    y=(double*)malloc(n*sizeof(double));
    a=b; k=1; r=1.0; z=dnmtcf(x,n);
    while (a>=eps)
      { for (i=0; i<=n-1; i++)
          y[i]=-a+2.0*a*mrnd1(&r)+x[i];
        z1=dnmtcf(y,n);
        k=k+1;
        if (z1>=z)
          { if (k>m) { k=1; a=a/2.0; }}
        else
          { k=1; 
            for (i=0; i<=n-1; i++) x[i]=y[i];
            z=z1;
            if (z<eps)  return;
          }
	}
    return;
}
#endif
/////////////////////////////////////////////////////////////
//用戶自編函數定義參考:
/*double dnmtcf(double x[],int n)
{
	double f,f1,f2,f3;
    n=n;
    f1=3.0*x[0]+x[1]+2.0*x[2]*x[2]-3.0;
    f2=-3.0*x[0]+5.0*x[1]*x[1]+2.0*x[0]*x[2]-1.0;
    f3=25.0*x[0]*x[1]+20.0*x[2]+12.0;
    f=sqrt(f1*f1+f2*f2+f3*f3);
    return(f);
}
double datknf(double x)
{
	double y;
    y=6.0-x*x;
    return(y);
}
double dpqrtf(double x)
{ 
	double y;
    y=x*x*(x-1.0)-1.0;
    return(y);
}
double ddhrtf(double x)
{ 
	double z;
    z=(((((x-5.0)*x+3.0)*x+1.0)*x-7.0)*x+7.0)*x-20.0;
    return(z);
}
void dnewtf(double x,double y[])
{ 
	y[0]=x*x*(x-1.0)-1.0;
    y[1]=3.0*x*x-2.0*x;
    return;
}
double dsnsef(double x[],double y[],int n)
{
	double z,f1,f2,f3,df1,df2,df3;
    n=n;
    f1=x[0]-5.0*x[1]*x[1]+7.0*x[2]*x[2]+12.0;
    f2=3.0*x[0]*x[1]+x[0]*x[2]-11.0*x[0];
    f3=2.0*x[1]*x[2]+40.0*x[0];
    z=f1*f1+f2*f2+f3*f3;
    df1=1.0; df2=3.0*x[1]+x[2]-11.0; df3=40.0;
    y[0]=2.0*(f1*df1+f2*df2+f3*df3);
    df1=10.0*x[1]; df2=3.0*x[0]; df3=2.0*x[2];
    y[1]=2.0*(f1*df1+f2*df2+f3*df3);
    df1=14.0*x[2]; df2=x[0]; df3=2.0*x[1];
    y[2]=2.0*(f1*df1+f2*df2+f3*df3);
    return(z);
}
void dnetnf(double x[],double y[],int n)
{ 
	y[0]=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]-1.0;
    y[1]=2.0*x[0]*x[0]+x[1]*x[1]-4.0*x[2];
    y[2]=3.0*x[0]*x[0]-4.0*x[1]+x[2]*x[2];
    n=n;
    return;
}
double dcmtcf(double x,double y)
{ 
	double u,v,z;
    u=x*x-y*y+x-y-2.0;
    v=2.0*x*y+x+y+2.0;
    z=sqrt(u*u+v*v);
    return(z);
}
double mrnd1(double* r)//From CH13 mrnd1.c
{ 
	int m;
    double s,u,v,p;
    s=65536.0; u=2053.0; v=13849.0;
    m=(int)(*r/s); *r=*r-m*s;
    *r=u*(*r)+v; m=(int)(*r/s);
    *r=*r-m*s; p=*r/s;
    return(p);
}
double dmtclf(double x)
{
	double y;
    y=exp(-x*x*x)-sin(x)/cos(x)+800.0;
    return(y);
}

void dnginf(int m,int n,double x[],double d[])
{ 
	m=m; n=n;
    d[0]=x[0]*x[0]+10.0*x[0]*x[1]+4.0*x[1]*x[1]+0.7401006;
    d[1]=x[0]*x[0]-3.0*x[0]*x[1]+2.0*x[1]*x[1]-1.0201228;
    return;
}

void dngins(int m,int n,double x[2],double* p)
{ m=m; n=n;
    p[0]=2.0*x[0]+7.0*x[1];
    p[1]=7.0*x[0]+6.0*x[1];
    p[2]=2.0*x[0]-2.0*x[1];
    p[3]=-2.0*x[0]+2.0*x[1];
    p[4]=1.0;
    p[5]=1.0;
    return;
}
*/
/////////////////////////////////////////////////////////////