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