今回は
(x-5)(x-4)(x-3)(x-2)(x-1)=0
を解く。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 5                             //方程式の最高次数
#define eps 1.0e-10             //精度決定定数

/*解の公式を用いて解を求める関数*/
void answer(double *a,double *b,double *c)
{
        double Re1,Re2,Im1,Im2,D;     
    
        D=*b**b-4.0**a**c;              //解の判別式                    
    
        if(D>0) {                           
                Re1=(-*b + sqrt(D))/(2.0**a);
                Re2=(-*b - sqrt(D))/(2.0**a);
                Im1=0.0;
                Im2=0.0;
        }
        else if(D==0) {
                Re1=-*b/(2.0**a);
                Im1=0.0;
        }
        else {                            
                Re1=Re2 =-*b/(2.0**a);
                Im1=sqrt(-D)/(2.0**a);
                Im2=-Im1;
        }

        if(D==0) {
                printf("%f + i %f\n",Re1,Im1);
        }
        else {
                printf("%f + i %f\n",Re1,Im1);
                printf("%f + i %f\n",Re2,Im2);
        }

}

/*ベアストウ法の主要な処理を行う関数*/
void bairstow(int n,double p[N+1])
{
        int i,max;
        double det,db,dc,q[N+1],u[N+1],w[N+1];
        double xxx,b,c;

        max=0;
        xxx=1.0;                //関数answerの引数*aに１を渡すため
        
        for(;;) {
                b=c=1.0;

                do {
                        max++;
                        if (max>1000000) {
                                printf("エラー\n");
                                exit(1);        //一定回数繰り返したらエラーとして終了とする
                        }

                        for(i=0;i<=n;i++) {
                                if(i==0) q[i]=p[i];
                        else if(i==1) q[i]=p[i]-b*q[i-1];
                                else q[i]=p[i]-b*q[i-1]-c*q[i-2];
                        }
                        for(i=0;i<=n;i++) {
                                if(i==0) u[i]=0.0;
                        else if(i==1) u[i]=-q[i-1]-b*u[i-1];
                                else u[i]=-q[i-1]-b*u[i-1]-c*u[i-2];
                        }
                        for(i=0;i<=n;i++) {
                                if(i==0) w[i]=0.0;
                        else if(i==1) w[i]=-b*w[i-1];
                        else w[i]=-q[i-2]-b*w[i-1]-c*w[i-2];
                        }

                        det=u[n-1]*w[n]-w[n-1]*(u[n]+q[n-1]);
                        db=1.0/det*(q[n]*w[n-1]-q[n-1]*w[n]);
                        dc=1.0/det*(-q[n]*u[n-1]+q[n-1]*(u[n]+q[n-1]));

                        b+=db;
                        c+=dc;
                } while ( db*db+dc*dc > eps*eps );

                for(i=0;i<=n-2;i++) p[i]=q[i];

                if(n>=3) {
                        answer(&xxx,&b,&c);
                        n-=2;
                }
                else if(n==2) {
                        answer(&p[0],&p[1],&p[2]);
                        break;
                }
                else if(n==1) {
                        printf("%f + i 0 \n",-p[1]/p[0]);
                        break;
                }
        }
}
        
int main()
{
        int n;
        double a[N+1]={1.0,-15.0,85.0,-225.0,274.0,-120.0};

        n=N;

        bairstow(n,a);

        return 0;
}

2.000000 + i 0.000000
1.000000 + i 0.000000
4.000000 + i 0.000000
3.000000 + i 0.000000
5.000000 + i 0