1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2002-2012, Industrial Light & Magic, a division of Lucas
4 // Digital Ltd. LLC
5 //
6 // All rights reserved.
7 //
8 // Redistribution and use in source and binary forms, with or without
9 // modification, are permitted provided that the following conditions are
10 // met:
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
16 // distribution.
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.
20 //
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.
32 //
33 ///////////////////////////////////////////////////////////////////////////
34
35
36
37 #ifndef INCLUDED_IMATHROOTS_H
38 #define INCLUDED_IMATHROOTS_H
39
40 //---------------------------------------------------------------------
41 //
42 // Functions to solve linear, quadratic or cubic equations
43 //
44 //---------------------------------------------------------------------
45
46 #include "ImathMath.h"
47 #include "ImathNamespace.h"
48 #include <complex>
49
50 IMATH_INTERNAL_NAMESPACE_HEADER_ENTER
51
52 //--------------------------------------------------------------------------
53 // Find the real solutions of a linear, quadratic or cubic equation:
54 //
55 // function equation solved
56 //
57 // solveLinear (a, b, x) a * x + b == 0
58 // solveQuadratic (a, b, c, x) a * x*x + b * x + c == 0
59 // solveNormalizedCubic (r, s, t, x) x*x*x + r * x*x + s * x + t == 0
60 // solveCubic (a, b, c, d, x) a * x*x*x + b * x*x + c * x + d == 0
61 //
62 // Return value:
63 //
64 // 3 three real solutions, stored in x[0], x[1] and x[2]
65 // 2 two real solutions, stored in x[0] and x[1]
66 // 1 one real solution, stored in x[1]
67 // 0 no real solutions
68 // -1 all real numbers are solutions
69 //
70 // Notes:
71 //
72 // * It is possible that an equation has real solutions, but that the
73 // solutions (or some intermediate result) are not representable.
74 // In this case, either some of the solutions returned are invalid
75 // (nan or infinity), or, if floating-point exceptions have been
76 // enabled with Iex::mathExcOn(), an Iex::MathExc exception is
77 // thrown.
78 //
79 // * Cubic equations are solved using Cardano's Formula; even though
80 // only real solutions are produced, some intermediate results are
81 // complex (std::complex<T>).
82 //
83 //--------------------------------------------------------------------------
84
85 template <class T> int solveLinear (T a, T b, T &x);
86 template <class T> int solveQuadratic (T a, T b, T c, T x[2]);
87 template <class T> int solveNormalizedCubic (T r, T s, T t, T x[3]);
88 template <class T> int solveCubic (T a, T b, T c, T d, T x[3]);
89
90
91 //---------------
92 // Implementation
93 //---------------
94
95 template <class T>
96 int
solveLinear(T a,T b,T & x)97 solveLinear (T a, T b, T &x)
98 {
99 if (a != 0)
100 {
101 x = -b / a;
102 return 1;
103 }
104 else if (b != 0)
105 {
106 return 0;
107 }
108 else
109 {
110 return -1;
111 }
112 }
113
114
115 template <class T>
116 int
solveQuadratic(T a,T b,T c,T x[2])117 solveQuadratic (T a, T b, T c, T x[2])
118 {
119 if (a == 0)
120 {
121 return solveLinear (b, c, x[0]);
122 }
123 else
124 {
125 T D = b * b - 4 * a * c;
126
127 if (D > 0)
128 {
129 T s = Math<T>::sqrt (D);
130 T q = -(b + (b > 0 ? 1 : -1) * s) / T(2);
131
132 x[0] = q / a;
133 x[1] = c / q;
134 return 2;
135 }
136 if (D == 0)
137 {
138 x[0] = -b / (2 * a);
139 return 1;
140 }
141 else
142 {
143 return 0;
144 }
145 }
146 }
147
148
149 template <class T>
150 int
solveNormalizedCubic(T r,T s,T t,T x[3])151 solveNormalizedCubic (T r, T s, T t, T x[3])
152 {
153 T p = (3 * s - r * r) / 3;
154 T q = 2 * r * r * r / 27 - r * s / 3 + t;
155 T p3 = p / 3;
156 T q2 = q / 2;
157 T D = p3 * p3 * p3 + q2 * q2;
158
159 if (D == 0 && p3 == 0)
160 {
161 x[0] = -r / 3;
162 x[1] = -r / 3;
163 x[2] = -r / 3;
164 return 1;
165 }
166
167 std::complex<T> u = std::pow (-q / 2 + std::sqrt (std::complex<T> (D)),
168 T (1) / T (3));
169
170 std::complex<T> v = -p / (T (3) * u);
171
172 const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
173 // for long double
174 std::complex<T> y0 (u + v);
175
176 std::complex<T> y1 (-(u + v) / T (2) +
177 (u - v) / T (2) * std::complex<T> (0, sqrt3));
178
179 std::complex<T> y2 (-(u + v) / T (2) -
180 (u - v) / T (2) * std::complex<T> (0, sqrt3));
181
182 if (D > 0)
183 {
184 x[0] = y0.real() - r / 3;
185 return 1;
186 }
187 else if (D == 0)
188 {
189 x[0] = y0.real() - r / 3;
190 x[1] = y1.real() - r / 3;
191 return 2;
192 }
193 else
194 {
195 x[0] = y0.real() - r / 3;
196 x[1] = y1.real() - r / 3;
197 x[2] = y2.real() - r / 3;
198 return 3;
199 }
200 }
201
202
203 template <class T>
204 int
solveCubic(T a,T b,T c,T d,T x[3])205 solveCubic (T a, T b, T c, T d, T x[3])
206 {
207 if (a == 0)
208 {
209 return solveQuadratic (b, c, d, x);
210 }
211 else
212 {
213 return solveNormalizedCubic (b / a, c / a, d / a, x);
214 }
215 }
216
217 IMATH_INTERNAL_NAMESPACE_HEADER_EXIT
218
219 #endif // INCLUDED_IMATHROOTS_H
220