|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
#include <ctime></ctime>
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
#include <cstdlib></cstdlib>
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
#include <cassert></cassert>
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
#include <complex></complex>
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
typedef double Real;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
typedef std::complex<real> Complex;</real>
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
const Real epsilon = 1e-6;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
const Real epsilonSqr = epsilon*epsilon;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
inline Real magSqr(Complex x)
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
{ return x.real()*x.real() + x.imag()*x.imag(); }
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
inline bool isZero(Real x)
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
{ return fabs(x) < epsilon; }
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
inline bool isZero(Complex x)
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
{ return magSqr(x) <= epsilonSqr; }
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
inline Complex verify2(Complex r, Real a, Real b, Real c) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
assert(isZero(a*r*r + b*r + c));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
return r;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
inline Complex verify3(Complex r, Real a, Real b, Real c, Real d) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
assert(isZero(a*r*r*r + b*r*r + c*r + d));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
return r;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
inline Complex verify4(Complex r, Real a, Real b, Real c, Real d, Real e) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
assert(isZero(a*r*r*r*r + b*r*r*r + c*r*r + d*r + e));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
return r;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
int solve2(Complex *roots, Real a, Real b, Real c) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if (isZero(a)) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if (isZero(b)) return 0;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[0] = verify2(-c/b, a, b, c);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
return 1;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real D = b*b - a*c*4;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real k = Real(0.5)/a;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex d = D < 0 ? Complex(0, sqrt(-D)) : Complex(sqrt(D));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[0] = verify2((-b - d)*k, a, b, c);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[1] = verify2((-b + d)*k, a, b, c);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
return 2;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
int solve3(Complex *roots, Real a, Real b, Real c, Real d) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if (isZero(a)) return solve2(roots, b, c, d);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// x = y - b/(3*a)
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// y*y*y + p*y + q = 0
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real p = (3*a*c - b*b)/(3*a*a);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real q = (27*a*a*d - 9*a*b*c + 2*b*b*b)/(27*a*a*a);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real Q = p*p*p/27 + q*q/4;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex Qs = Q < 0 ? Complex(0, sqrt(-Q)) : Complex(sqrt(Q));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex A = pow(-q/2 + Qs, Real(1)/3);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex B = pow(-q/2 - Qs, Real(1)/3);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// choose complimentary B for A (A*B must be equal -p/3)
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex rot(Real(-0.5), sqrt(Real(3))/2);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if (!isZero(A*B + p/3)) B *= rot;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if (!isZero(A*B + p/3)) B *= rot;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex Y = (A - B)*Complex(0, sqrt(Real(3)));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex y0 = A + B;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex y1 = (-y0 - Y)/Real(2);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex y2 = (-y0 + Y)/Real(2);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real dd = b/(3*a);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[0] = verify3(y0 - dd, a, b, c, d);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[1] = verify3(y1 - dd, a, b, c, d);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[2] = verify3(y2 - dd, a, b, c, d);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
return 3;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
int solve4(Complex *roots, Real a, Real b, Real c, Real d, Real e) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if (isZero(a)) return solve3(roots, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
//printf("%g*x^4 + %g*x^3 + %g*x^2 + %g*x + %g = 0\n", a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// x = y - b/(4*a)
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// y^4 + p*y^2 + q*y + r = 0
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real dd = b/(4*a);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real p = (8*a*c - 3*b*b)/(8*a*a);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real q = (8*a*a*d - 4*a*b*c + b*b*b)/(8*a*a*a);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real r = (256*a*a*a*e - 64*a*a*b*d + 16*a*b*b*c - 3*b*b*b*b)/(256*a*a*a*a);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
//printf(", y^4 + %g*y^2 + %g*y + %g = 0", p, q, r);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if (isZero(q)) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// biquadratic equation
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// y^4 + p*y^2 + r = 0
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex y[2];
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
solve2(y, 1, p, r);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
y[0] = sqrt(y[0]);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
y[1] = sqrt(y[1]);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[0] = verify4(-y[0] - dd, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[1] = verify4( y[0] - dd, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[2] = verify4(-y[1] - dd, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[3] = verify4( y[1] - dd, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
return 4;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// solve cubic equation
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// z*z*z + (p/2)*z*z + ((p*p - 4*r)/16)*z - q*q/64 = 0
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real pp = p/2;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real qq = (p*p - 4*r)/16;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real rr = -q*q/64;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
//printf("z^3 + %g*z^2 + %g*z + %g = 0\n", pp, qq, rr);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex z[3];
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
solve3(z, 1, pp, qq, rr);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
z[0] = sqrt(z[0]);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
z[1] = sqrt(z[1]);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
z[2] = sqrt(z[2]);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// we need to find signs combination where following is valid:
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// (+-z0)*(+-z1)*(+-z2) = -q/8
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex zzz = z[0]*z[1]*z[2];
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
//printf(
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// "sqrt(z): (%g, %g), (%g, %g), (%g, %g), zzz: (%g, %g)\n",
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// z[0].real(), z[0].imag(),
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// z[1].real(), z[1].imag(),
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// z[2].real(), z[2].imag(),
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
// zzz.real(), zzz.imag() );
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
assert(isZero(zzz.imag()));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
assert(isZero(fabs(zzz.real()) - fabs(q/8)));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if ((zzz.real() > 0) == (q > 0))
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
z[0] = -z[0];
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
assert(isZero(z[0]*z[1]*z[2] + q/8));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[0] = verify4( z[0] - z[1] - z[2] - dd, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[1] = verify4(-z[0] + z[1] - z[2] - dd, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[2] = verify4(-z[0] - z[1] + z[2] - dd, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
roots[3] = verify4( z[0] + z[1] + z[2] - dd, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
return 4;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
int main() {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
srand(time(NULL));
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
bool failed = false;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
for(int i = 0; i < 1000000; ++i) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real a = rand()%20 - 10;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real b = rand()%20 - 10;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real c = rand()%20 - 10;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real d = rand()%20 - 10;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Real e = rand()%20 - 10;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
printf("%6d: %g*x^4 + %g*x^3 + %g*x^2 + %g*x + %g = 0", i, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
int vcnt = a ? 4
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
: b ? 3
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
: c ? 2
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
: d ? 1 : 0;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex roots[4];
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
int cnt = solve4(roots, a, b, c, d, e);
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
bool f = cnt != vcnt;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
for(int i = 0; i < cnt; ++i) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
const Complex &x = roots[i];
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
printf(", (%g, %g)", x.real(), x.imag());
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
Complex r = a*x*x*x*x + b*x*x*x + c*x*x + d*x + e;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if (!isZero(r))
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
{ printf(" ![%g, %g]", r.real(), r.imag()); f = true; }
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
printf(", %s\n", f ? "FAILED" : "passed");
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
if (f) {
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
printf("^^^^^^^^^^\n");
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
failed = true;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
break;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
printf("\n%s\n\n", failed ? "FAILED" : "success");
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
return failed;
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
}
|
|
![](https://seccdn.libravatar.org/avatar/2e5dd0bee1e7e619066117de357c8458fc7e847f4345b0cb8a7a5413aa2a45a0?s=16&d=retro) |
62033a |
|