多项式计算的Horner 方法

Horner 算法是以英国数学家 William George Horner 命名的一种多项式求值的快速算法，实际上，这种快速算法在他之前就已经被Paolo Ruffini使用过了。而中国数学家秦九韶提出这种算法要比William George Horner 早600多年。

P(x) 是一个多项式：

我们希望计算x取某个特殊值x0时多项式的值p(x0).

构造一个序列：

那么这个序列b0的值就是多项式的值了。

用程序实现如下：

double horner(double p[], int n, double x)
{
    double sum;
    sum = p[--n];
    while ( n > 0 )
    {
        sum = p[--n] + sum * x;
    }
    return sum;
}

我经常要设计FIR、IIR 滤波器，设计完滤波器后都要验证一下频率特性是否正确。这时就需要计算实系数多项式在x取复数值时的结果。因此就有了下面的代码：

 double _Complex horner_C(double p[], int n, double _Complex x)
 {
    double _Complex sum;
    sum = p[--n];
    while ( n > 0 )
    {
        sum = p[--n] + sum * x;
    }
    return sum;
 }

这个代码用到了C99中对复数支持的新特性，需要编译器支持C99标准。

最后给一个测试代码，其中多项式p构成了一个FIR低通滤波器，这个低通滤波器是我用scilab 中的ffilt 函数设计的。

滤波器构造的代码如下：

ffilt('lp',20 , 0.2 , 0.5);

下面是计算这个滤波器的频响特性的代码。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
double horner(double p[], int n, double x);
double _Complex horner_C(double p[], int n, double _Complex x);
int main()
{
    int i;
    double p[] = {
    -0.0196945,    -0.0356154,
    1.559E-17,    0.0465740,
    0.0340178,    -0.0415773,
    -0.0864945,    1.559E-17,
    0.2018205,    0.3741957,
    0.3741957,    0.2018205,
    1.559E-17,    -0.0864945,
    -0.0415773,    0.0340178,
    0.0465740,    1.559E-17,
    -0.0356154,    -0.0196945
    };
    double x[201], mag[201], pha[201];
    double _Complex xx;
    for(i = 0; i <= 200; i++)
    {
        x[i] = 0.5 * i / 200.0;
        xx = cexp(-2 * M_PI * x[i] * I);
        mag[i] = cabs(horner_C(p, 20, xx));
        pha[i] = carg(horner_C(p, 20, xx));
        printf("%f, %f, %f\n", x[i] , mag[i], pha[i]);
    }

    return 0;
}

将输出结果画出了可以得到如下的幅频特性：

scilab 计算出的幅频特性如下图：

可以看出两幅图完全相同，说明我们的代码计算结果是正确的。