HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ImathRoots.h
Go to the documentation of this file.
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 
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
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
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
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
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 
218 
219 #endif // INCLUDED_IMATHROOTS_H
int solveNormalizedCubic(T r, T s, T t, T x[3])
Definition: ImathRoots.h:151
#define IMATH_INTERNAL_NAMESPACE_HEADER_EXIT
#define IMATH_INTERNAL_NAMESPACE_HEADER_ENTER
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER int solveLinear(T a, T b, T &x)
Definition: ImathRoots.h:97
const GLdouble * v
Definition: glcorearb.h:836
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
SYS_API double pow(double x, double y)
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
int solveCubic(T a, T b, T c, T d, T x[3])
Definition: ImathRoots.h:205
GLint GLenum GLint x
Definition: glcorearb.h:408
GLboolean r
Definition: glcorearb.h:1221
static T sqrt(T x)
Definition: ImathMath.h:115
int solveQuadratic(T a, T b, T c, T x[2])
Definition: ImathRoots.h:117