多项式与连分式函数的计算

#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "ComplexNumber.h"
#include "Multinomial.h"

//复系数多项式求值
void ncply(double ar[], double ai[], int n, double x, double y, double *u, double *v)
{
 int i;
 double p,q,s,t;
 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 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;
}

//复系数多项式相乘
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;
 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 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 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;
 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 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);
}

//函数连分式的计算
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);
}

//二维多项式求值
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);

}

                                 ----根据《C语言常用算法程序集》整理