/************************************************
 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 "CH14.h"  多項式與連分式函數(shù)的計算*****/
#ifndef CH14_H_
#define CH14_H_

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

#include "ch15.h"
//*******************************************************************

double nplyv(double a[],int n,double x);//一維多項式求值
void nplys(double a[],int n,double x[],int m,double p[]);//一維多項式多組求值
double nbply(double a[],int m,int n,double x,double y);//二維多項式求值
void ncply(double ar[],double ai[],int n,double x,double y,double* u,double* v);//復系數(shù)多項式求值
void npmul(double p[],int m,double q[],int n,double s[],int k);//多項式相乘
void npdiv(double p[],int m,double q[],int n,double s[],int k,double r[],int l);//多項式相除
void ncmul(double pr[],double pi[],int m,double qr[],
		   double qi[],int n,double sr[],double si[],int k);//復系數(shù)多項式相乘
void ncdiv(double pr[],double pi[],int m,double qr[],double qi[],int n,double sr[],
		   double si[],int k,double rr[],double ri[],int l);//復系數(shù)多項式相除
double nfpqv(double x[],double b[],int n,double t);//函數(shù)連分式的的計算

//*******************************************************************
double nplyv(double a[],int n,double x)
{ 
	int i;
    double u;
    u=a[n-1];
    for (i=n-2; i>=0; i--)
      u=u*x+a[i];
    return(u);
}
/////////////////////////////////////////////////////////////
void nplys(double a[],int n,double x[],int m,double p[])
{ 
	int i,j,mm,nn,ll,t,s,kk,k;
    double *b,y,z;
    b=(double*)malloc(2*n*sizeof(double));
    y=a[n-1];
    for (i=0; i<=n-1; i++) b[i]=a[i]/y;
    k=log(n-0.5)/log(2.0)+1;
    nn=1;
    for (i=0; i<=k-1; i++) nn=2*nn;
    for (i=n; i<nn-1; i++) b[i]=0.0;
    b[nn-1]=1.0;
    t=nn; s=1;
    for (i=1; i<=k-1; i++)
      { t=t/2; mm=-t;
        for (j=1; j<=s; j++)
          { mm=mm+t+t; b[mm-1]=b[mm-1]-1.0; 
            for (kk=2; kk<=t; kk++)
              b[mm-kk]=b[mm-kk]-b[mm-1]*b[mm+t-kk];
          }
        s=s+s;
      }
    for (kk=1; kk<=m; kk++)
      { for (i=0; i<=(nn-2)/2; i++)
           a[i]=x[kk-1]+b[2*i];
        mm=1; z=x[kk-1];
        for (i=1; i<=k-1; i++)
          { mm=mm+mm; ll=mm+mm; z=z*z;
            for (j=0; j<=nn-1; j=j+ll)
              a[j/2]=a[j/2]+a[(j+mm)/2]*(z+b[j+mm-1]);
          }
        z=z*z/x[kk-1];
        if (nn!=n) a[0]=a[0]-z;
        p[kk-1]=a[0]*y;
      }
    return;
}
/////////////////////////////////////////////////////////////
double nbply(double a[],int m,int n,double x,double y)
{ 
	int i,j;
    double u,s,xx;
    u=0.0; xx=1.0;
    for (i=0; i<=m-1; i++)
      { s=a[i*n+n-1]*xx;
        for (j=n-2; j>=0; j--)
          s=s*y+a[i*n+j]*xx;
        u=u+s; xx=xx*x;
      }
    return(u);
}
/////////////////////////////////////////////////////////////
void ncply(double ar[],double ai[],int n,double x,double y,double* u,double* v)
{ 
	int i;
    double p,q,s,t;
    //extern void ocmul();
    s=ar[n-1]; t=ai[n-1];
    for (i=n-2; i>=0; i--)
      { ocmul(s,t,x,y,&p,&q);
        s=p+ar[i]; t=q+ai[i];
      }
    *u=s; *v=t;
    return;
}
/////////////////////////////////////////////////////////////
void npmul(double p[],int m,double q[],int n,double s[],int k)
{ 
	int i,j;
    for (i=0; i<=k-1; i++) s[i]=0.0;
    for (i=0; i<=m-1; i++)
    for (j=0; j<=n-1; j++)
      s[i+j]=s[i+j]+p[i]*q[j];
    return;
}
/////////////////////////////////////////////////////////////
void npdiv(double p[],int m,double q[],int n,double s[],int k,double r[],int l)
{ 
	int i,j,mm,ll;
    for (i=0; i<=k-1; i++) s[i]=0.0;
    if (q[n-1]+1.0==1.0) return;
    ll=m-1;
    for (i=k; i>=1; i--)
      { s[i-1]=p[ll]/q[n-1];
        mm=ll;
        for (j=1; j<=n-1; j++)
          { p[mm-1]=p[mm-1]-s[i-1]*q[n-j-1];
            mm=mm-1;
          }
        ll=ll-1;
      }
    for (i=0; i<=l-1; i++) r[i]=p[i];
    return;
}
/////////////////////////////////////////////////////////////
void ncmul(double pr[],double pi[],int m,double qr[],double qi[],int n,double sr[],double si[],int k)
{ 
	int i,j;
    double a,b,c,d,u,v;
    //extern void ocmul();
    for (i=0; i<=k-1; i++)
      { sr[i]=0.0; si[i]=0.0; }
    for (i=0; i<=m-1; i++)
    for (j=0; j<=n-1; j++)
      { a=pr[i]; b=pi[i]; c=qr[j]; d=qi[j];
        ocmul(a,b,c,d,&u,&v);
        sr[i+j]=sr[i+j]+u;
        si[i+j]=si[i+j]+v;
      }
    return;
}
/////////////////////////////////////////////////////////////
void ncdiv(double pr[],double pi[],int m,double qr[],double qi[],
		   int n,double sr[],double si[],int k,double rr[],double ri[],int l)
{ 
	int i,j,mm,ll;
    double a,b,c,d,u,v;
    //extern void ocmul();
    //extern void ocdiv();
    for (i=0; i<=k-1; i++) 
      { sr[i]=0.0; si[i]=0.0; }
    d=qr[n-1]*qr[n-1]+qi[n-1]*qi[n-1];
    if (d+1.0==1.0) return;
    ll=m-1;
    for (i=k; i>=1; i--)
      { a=pr[ll]; b=pi[ll]; c=qr[n-1]; d=qi[n-1];
        ocdiv(a,b,c,d,&u,&v);
        sr[i-1]=u; si[i-1]=v;
        mm=ll;
        for (j=1; j<=n-1; j++)
          { a=sr[i-1]; b=si[i-1]; 
            c=qr[n-j-1]; d=qi[n-j-1];
            ocmul(a,b,c,d,&u,&v);
            pr[mm-1]=pr[mm-1]-u;
            pi[mm-1]=pi[mm-1]-v;
            mm=mm-1;
          }
        ll=ll-1;
      }
    for (i=0; i<=l-1; i++) 
      { rr[i]=pr[i]; ri[i]=pi[i]; }
    return;
}
/////////////////////////////////////////////////////////////
double nfpqv(double x[],double b[],int n,double t)
{
	int k;
    double u;
    u=b[n-1];
    for (k=n-2; k>=0; k--)
      { if (fabs(u)+1.0==1.0)
           u=1.0e+35*(t-x[k])/fabs(t-x[k]);
        else
           u=b[k]+(t-x[k])/u;
      }
    return(u);
}
/////////////////////////////////////////////////////////////
#endif