48{
49 double m, c1, c0;
50
51
52
53
54
55 std::complex<double> de = A[0][1] * A[1][2];
59 m = real(A[0][0]) + real(A[1][1]) + real(A[2][2]);
60 c1 = (real(A[0][0])*real(A[1][1])
61 + real(A[0][0])*real(A[2][2])
62 + real(A[1][1])*real(A[2][2]))
63 - (dd + ee + ff);
64 c0 = real(A[2][2])*dd + real(A[0][0])*ee + real(A[1][1])*ff
65 - real(A[0][0])*real(A[1][1])*real(A[2][2])
66 - 2.0 * (real(A[0][2])*real(de) + imag(A[0][2])*imag(de));
67
68
69 double p, sqrt_p, q, c, s, phi;
71 q = m*(p - (3.0/2.0)*c1) - (27.0/2.0)*c0;
72 sqrt_p = sqrt(fabs(p));
73
74 phi = 27.0 * ( 0.25*
SQR(c1)*(p - c1) + c0*(q + 27.0/4.0*c0));
75 phi = (1.0/3.0) * atan2(sqrt(fabs(phi)), q);
76
77 c = sqrt_p*cos(phi);
78 s = (1.0/
M_SQRT3)*sqrt_p*sin(phi);
79
80 w[1] = (1.0/3.0)*(m - c);
81 w[2] = w[1] + s;
82 w[0] = w[1] + c;
83 w[1] -= s;
84
85 return 0;
86}