1 ///////////////////////////////////////////////////////////////////////////
3 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
6 // All rights reserved.
8 // Redistribution and use in source and binary forms, with or without
9 // modification, are permitted provided that the following conditions are
11 // * Redistributions of source code must retain the above copyright
12 // notice, this list of conditions and the following disclaimer.
13 // * Redistributions in binary form must reproduce the above
14 // copyright notice, this list of conditions and the following disclaimer
15 // in the documentation and/or other materials provided with the
17 // * Neither the name of Industrial Light & Magic nor the names of
18 // its contributors may be used to endorse or promote products derived
19 // from this software without specific prior written permission.
21 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 ///////////////////////////////////////////////////////////////////////////
37 #ifndef INCLUDED_IMATHROOTS_H
38 #define INCLUDED_IMATHROOTS_H
40 //---------------------------------------------------------------------
42 // Functions to solve linear, quadratic or cubic equations
44 //---------------------------------------------------------------------
50 //--------------------------------------------------------------------------
51 // Find the real solutions of a linear, quadratic or cubic equation:
53 // function equation solved
55 // solveLinear (a, b, x) a * x + b == 0
56 // solveQuadratic (a, b, c, x) a * x*x + b * x + c == 0
57 // solveNormalizedCubic (r, s, t, x) x*x*x + r * x*x + s * x + t == 0
58 // solveCubic (a, b, c, d, x) a * x*x*x + b * x*x + c * x + d == 0
62 // 3 three real solutions, stored in x[0], x[1] and x[2]
63 // 2 two real solutions, stored in x[0] and x[1]
64 // 1 one real solution, stored in x[1]
65 // 0 no real solutions
66 // -1 all real numbers are solutions
70 // * It is possible that an equation has real solutions, but that the
71 // solutions (or some intermediate result) are not representable.
72 // In this case, either some of the solutions returned are invalid
73 // (nan or infinity), or, if floating-point exceptions have been
74 // enabled with Iex::mathExcOn(), an Iex::MathExc exception is
77 // * Cubic equations are solved using Cardano's Formula; even though
78 // only real solutions are produced, some intermediate results are
79 // complex (std::complex<T>).
81 //--------------------------------------------------------------------------
83 template <class T> int solveLinear (T a, T b, T &x);
84 template <class T> int solveQuadratic (T a, T b, T c, T x[2]);
85 template <class T> int solveNormalizedCubic (T r, T s, T t, T x[3]);
86 template <class T> int solveCubic (T a, T b, T c, T d, T x[3]);
95 solveLinear (T a, T b, T &x)
115 solveQuadratic (T a, T b, T c, T x[2])
119 return solveLinear (b, c, x[0]);
123 T D = b * b - 4 * a * c;
129 x[0] = (-b + s) / (2 * a);
130 x[1] = (-b - s) / (2 * a);
148 solveNormalizedCubic (T r, T s, T t, T x[3])
150 T p = (3 * s - r * r) / 3;
151 T q = 2 * r * r * r / 27 - r * s / 3 + t;
154 T D = p3 * p3 * p3 + q2 * q2;
156 if (D == 0 && p3 == 0)
164 std::complex<T> u = std::pow (-q / 2 + std::sqrt (std::complex<T> (D)),
167 std::complex<T> v = -p / (T (3) * u);
169 const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
171 std::complex<T> y0 (u + v);
173 std::complex<T> y1 (-(u + v) / T (2) +
174 (u - v) / T (2) * std::complex<T> (0, sqrt3));
176 std::complex<T> y2 (-(u + v) / T (2) -
177 (u - v) / T (2) * std::complex<T> (0, sqrt3));
181 x[0] = y0.real() - r / 3;
186 x[0] = y0.real() - r / 3;
187 x[1] = y1.real() - r / 3;
192 x[0] = y0.real() - r / 3;
193 x[1] = y1.real() - r / 3;
194 x[2] = y2.real() - r / 3;
202 solveCubic (T a, T b, T c, T d, T x[3])
206 return solveQuadratic (b, c, d, x);
210 return solveNormalizedCubic (b / a, c / a, d / a, x);