这里可以将这两个函数分别对double数组和complex数组实现
double_filepud
void double_filepud(double *m, int length)
{
int i, j;
double temp;
for (i = 0; i < length / 2; i++)
{
// 交换当前元素和对应的对称元素
temp = m[i];
m[i] = m[length - 1 - i];
m[length - 1 - i] = temp;
}
}
由于filepud函数比较简单,就不罗列complex数组的filepud实现了
double_poly
void double_poly(double *p, int np, double *d)
{
double *temp = (double *)malloc(sizeof(double) * np * np);
for (size_t i = 0; i < np; i++)
{
for (size_t j = 0; j < np; j++)
{
temp[i * np + j] = p[i * np + j];
}
}
// eig(A)
double *u_real = (double *)malloc(sizeof(double) * np);
double *v_imag = (double *)malloc(sizeof(double) * np);
int jt = 60;
hhbg(temp, np);
hhqr(temp, np, EPS, jt, u_real, v_imag);
Complex *dt = (Complex *)malloc(sizeof(Complex) * (np + 1));
Complex *dtt = (Complex *)malloc(sizeof(Complex) * (np + 1));
dt[0].real = 1;
dt[0].imaginary = 0;
for (size_t i = 1; i < np + 1; i++)
{
dt[i].real = 0;
dt[i].imaginary = 0;
}
for (size_t i = 0; i < np; i++)
{
for (size_t j = 1; j <= i + 1; j++)
{
dtt[j] = complex_subtract(dt[j], complex_multiply(create(u_real[i], v_imag[i]), dt[j - 1]));
}
for (size_t j = 1; j <= i + 1; j++)
{
dt[j].real = dtt[j].real;
dt[j].imaginary = dtt[j].imaginary;
}
}
for (size_t i = 0; i < np + 1; i++)
{
d[i] = dt[i].real;
}
free(temp);
free(u_real);
free(v_imag);
free(dt);
free(dtt);
}
这个函数针对矩阵,先转换为上海森伯格矩阵,再使用QR法来求解
complex_poly
void complex_poly(Complex *p, int np, Complex *d)
{
Complex *c = (Complex *)malloc((np + 1) * sizeof(Complex));
c[0].real = 1.0;
c[0].imaginary = 0.0;
d[0] = c[0];
for (int i = 1; i < np + 1; i++)
{
c[i].real = 0.0;
c[i].imaginary = 0.0;
d[i]=c[i];
}
Complex temp;
for (int i = 0; i < np; i++)
{
for (int j = 1; j <= i + 1; j++)
{
temp = complex_multiply(p[i], d[j - 1]);
c[j].real = d[j].real - temp.real;
c[j].imaginary = d[j].imaginary - temp.imaginary;
}
for (int j = 1; j <= i + 1; j++)
{
d[j].real = c[j].real;
d[j].imaginary = c[j].imaginary;
}
}
free(c);
}
在matlab中向量也可以使用poly,这里再加上针对double向量的ploy函数实现
vector_poly
void vector_poly(double *p, int np, double *d)
{
// 初始化多项式系数
double *dt = (double *)malloc(sizeof(double) * (np + 1));
double *dtt = (double *)malloc(sizeof(double) * (np + 1));
for (int i = 0; i <= np; ++i) {
dt[i] = 0.0;
}
dt[0] = 1.0;
for (size_t i = 0; i < np; i++)
{
for (size_t j = 1; j <= i + 1; j++)
{
dtt[j] = dt[j]-p[i]*dt[j - 1];
}
for (size_t j = 1; j <= i + 1; j++)
{
dt[j] = dtt[j];
}
}
for (size_t i = 0; i < np + 1; i++)
{
d[i] = dt[i];
}
free(dt);
free(dtt);
}
测试
经过与matlab结果进行对比,基本无误,可以放心食用