3x3 実対称行列固有値

を計算する超手抜きコードを30分で書いた。
Det|M-xI|=0の三次方程式をニュートン法で解いてる。初期値はまんなかの変曲点 -c2/3 と、微分０になる二箇所 -c2/3±sqrt(c2*c2-3*c1)/3の外側。解がないとか収束判定とかまったく考えてない。

static double newton3(double x0, double c0, double c1, double c2)
{
    double x = x0;
    int i;

    for(i=0;i<1000;i++)
    {
	double f = x*x*x + c2*x*x+c1*x +c0;
	double df = 3*x*x + 2*c2*x + c1;

	x -= f/df;
    }
    return x;
}

extern void mat_eig3(double m[3][3], double v[3])
{
    double m11=m[0][0];
    double m22=m[1][1];
    double m33=m[2][2];
    double m12=m[0][1];
    double m23=m[1][2];
    double m31=m[2][0];

    double c0 =   -m11*m22*m33 
	        -2*m12*m23*m31 
		  +m11*m23*m23
		  +m22*m31*m31
		  +m33*m12*m12;

    double c1 = m22*m33+m33*m11+m11*m22 - (m23*m23 + m12*m12 + m31*m31);

    double c2 = -(m11+m22+m33);

    v[1] = newton3( -c2/3.0, c0, c1, c2);
    v[0] = newton3( -c2/3.0 - 2*sqrt(c2*c2-3*c1)/3, c0, c1, c2);
    v[2] = newton3( -c2/3.0 + 2*sqrt(c2*c2-3*c1)/3, c0, c1, c2);
}