多项式求值 python/c 程序
#!coding: utf-8
#多项式求值，多组求值法
#F(x)首一化后变为T(x)，T(x)是2^k项的多项式，分解为(x^(k-1)+b)*q(x)+r(x)，其中q(x)和r(x)都是2^(k-1)项的多项式
#根据 常用算法程序集（徐士良编，ISBN：730209697) c源代码/ch1/1plys.c改写
#基本借鉴了1plys.c程序的两两合并思想，但去掉了系数首一化，因此，系数预处理只是对
#非2的次方的多项式系数扩充到2的次方，没有做其他系数转换
#求k的过程更加易懂，见log_2程序
#求值时两两合并用的是递归程序，体现出次算法的递归性质，代码更简单明白
#实验表明，计算结果精确无误
#以下是c算法的结果：
#[ning@58 ~]$ ./desktop/notebook/material/20_10_2008/a.out
#x(0)= 0.90 p(0)=-1.8562268e+01
#x(1)=-0.90 p(1)=-2.6715368e+01
#x(2)= 1.10 p(2)=-1.9556128e+01
#x(3)=-1.10 p(3)=-2.1513028e+01
#x(4)= 1.30 p(4)=-2.0875732e+01
#x(5)=-1.30 p(5)=-6.3404320e+00
#以下是本程序的结果：
#[ning@58 python_math_lab]$ python polynomial_v.py
#[-18.562268]
#[-26.715367999999998]
#[-19.556128000000001]
#[-21.513027999999998]
#[-20.875731999999999]
#[-6.3404319999999927]
def coefficient_pretreatment(a):
(k,nn)=log_2(len(a))
a+=(nn-len(a))*[0]
return (a,k,nn)
def log_2(n):
k=nn=1
nn=1
while nn n:
nn=1
k+=1
return (k,nn)
def recursive_poly_v(a,k,nn,x):
if k>=1:
a=[a+a[i+1]*x for i in xrange(0,nn,2)]
k -= 1
nn >>= 1
x=x*x
recursive_poly_v(a,k,nn,x)
else:
print a
a=[-20.0,7.0,-7.0,1.0,3.0,-5.0,2.0]
(a,k,nn)=coefficient_pretreatment(a)
x=[0.9,-0.9,1.1,-1.1,1.3,-1.3]
for i in x:
recursive_poly_v(a,k,nn,i)
以下是1plys.c，注释是我加的/*
* 注意：链接math.h库，用gcc编译时，要手动加上-lm参数
用系数与处理法对多项式求值.
其中n=2^k(k>=1)
首先将多项式p(x)变为首一多项式T(x),然后对系数进行预处理
预处理过程:
将T(x)分解为
T(x)=(x^j+b)*q(x)+r(x)
其中j=2^(k-1)
b=T[2^(k-1)-1]-1
q(x)的系数是T(x)前一半,
r(x)的系数是T(x)后一半减去b*q(x)
这样,使得 q(x)与r(x)均为2^(k-1)-1次的首一多项式
然后对q(x)和r(x)进行系数预处理
继续上述过程,直到q'(x)成为1,而r'(x)等于0为止
double a[n] 存放n-1次多项式的n个系数. 返回时其值将被改变
int n 多项式的项数
double x[m] 存放给定的m个自变量值
int m 给定自变量个数
double p[m] 返回时存放与给定m个自变量值对应的多项式值
void plys() 过程
*/
#include math.h>
#include stdlib.h>
int plys(a,n,x,m,p)
int n,m;
double a[],x[],p[];
{ int i,j,mm,nn,ll,t,s,kk,k;
double *b,y,z;
b=malloc(2*n*sizeof(double)); /*b的空间是a的2倍*/
y=a[n-1];
for (i=0; i=n-1; i++) b=a/y; /*多项式首一化*/
k=log(n-0.5)/log(2.0)+1; /*得到k值，至于为何要先减去0.5，后再加1，请看下面用python做的实验*/
nn=1;
for (i=0; i=k-1; i++) nn=2*nn; /*得到2^k, nn是扩充后的项数*/
for (i=n; inn-1; i++) b=0.0; /*多项式项数不是2^k, 则用0补足*/
b[nn-1]=1.0; /*????????????*/
t=nn; s=1; /*t是每段长度, s是奇数段的个数*/
for (i=1; i=k-1; i++)
{ t=t/2; mm=-t; /*多项式被切成2^i截,mm始终为第奇数段末尾,t是每段的长度*/
for (j=1; j=s; j++) /*现在处理第j个奇数段*/
{ mm=mm+t+t; b[mm-1]=b[mm-1]-1.0; /*mm跳到下一个奇数段, b[mm-1]就是传说中的b*/
for (kk=2; kk=t; kk++) /*b[后半部分] -= b*b[前半部分]*/
b[mm-kk]=b[mm-kk]-b[mm-1]*b[mm+t-kk];
}
s=s+s; /*s变成2, 因为多项式被划成两截*/
}
/*
-----------------------------------------------------------------------
以上，系数预处理完成，处理结果放在b数组中
以下，对m个自变量求多项式值，m是自变量个数，nn是2的k次方
*/
for (kk=1; kk=m; kk++)
{ for (i=0; i=(nn-2)/2; i++)/*做第一次迭代,并把结果保存在a中; 此时是1次多项式*/
a=x[kk-1]+b[2*i]; /*a奇数项是所谓的r(x),偶数项是所谓的q(x)*/
mm=1; z=x[kk-1];
for (i=1; i=k-1; i++) /*合并k-1次,把2的k-1次项合并成一个，i控制合并的跨度，一开始是4，体现在下面的ll上*/
{ mm=mm+mm; ll=mm+mm; z=z*z; /*mm一开始是2，ll一开始是4；开始时在b数组上4个4个前进，这时x的次方是2，下次循环要加倍*/
for (j=0; j=nn-1; j=j+ll) /*从左到右进行i合并*/
a[j/2]=a[j/2]+a[(j+mm)/2]*(z+b[j+mm-1]);/*b[j+mm-1]是对应的b; q = q + r*(x^mm+b)r, 结果放置在前面操作数，这样最后结果在a[0]中*/
}
z=z*z/x[kk-1]; /*z的nn－1次*/
if (nn!=n) a[0]=a[0]-z; /*销去b[nn-1]=1.0带来的影响*/
p[kk-1]=a[0]*y; /*反首一化*/
}
return 0;
}
/*
在注释中有如此多形如'x的xx', 可以感悟出面向对象的思想是从何产生,
面向对象使程序更具可读性
*/
/*
python 实验：
>>> import math
>>> def t(x):
... return math.log(x-0.5)/math.log(2.0)+1
...
>>> int(t(4))
2
>>> int(t(5))
3
>>> int(t(3))
2
>>> def t2(x):
... return math.log(x)/math.log(2.0)
...
>>> int(t2(4))
2
>>> int(t2(3))
1
>>> int(t2(5))
2
>>>
*/
