matlab flipud poly函数 C实现

• filpud(数组从上向下翻转)

• poly(具有指定根的多项式或特征多项式)

这里可以将这两个函数分别对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结果进行对比，基本无误，可以放心食用