source: NonGTP/OpenEXR/include/Imath/ImathRoots.h @ 855

Revision 855, 5.5 KB checked in by igarcia, 19 years ago (diff)
RevLine 
[855]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 <complex>
47
48namespace Imath {
49
50//--------------------------------------------------------------------------
51// Find the real solutions of a linear, quadratic or cubic equation:
52//
53//      function                                   equation solved
54//
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
59//
60// Return value:
61//
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
67//
68// Notes:
69//
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
75//      thrown.
76//
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>).
80//
81//--------------------------------------------------------------------------
82
83template <class T> int  solveLinear (T a, T b, T &x);
84template <class T> int  solveQuadratic (T a, T b, T c, T x[2]);
85template <class T> int  solveNormalizedCubic (T r, T s, T t, T x[3]);
86template <class T> int  solveCubic (T a, T b, T c, T d, T x[3]);
87
88
89//---------------
90// Implementation
91//---------------
92
93template <class T>
94int
95solveLinear (T a, T b, T &x)
96{
97    if (a != 0)
98    {
99        x = -b / a;
100        return 1;
101    }
102    else if (b != 0)
103    {
104        return 0;
105    }
106    else
107    {
108        return -1;
109    }
110}
111
112
113template <class T>
114int
115solveQuadratic (T a, T b, T c, T x[2])
116{
117    if (a == 0)
118    {
119        return solveLinear (b, c, x[0]);
120    }
121    else
122    {
123        T D = b * b - 4 * a * c;
124
125        if (D > 0)
126        {
127            T s = sqrt (D);
128
129            x[0] = (-b + s) / (2 * a);
130            x[1] = (-b - s) / (2 * a);
131            return 2;
132        }
133        if (D == 0)
134        {
135            x[0] = -b / (2 * a);
136            return 1;
137        }
138        else
139        {
140            return 0;
141        }
142    }
143}
144
145
146template <class T>
147int
148solveNormalizedCubic (T r, T s, T t, T x[3])
149{
150    T p  = (3 * s - r * r) / 3;
151    T q  = 2 * r * r * r / 27 - r * s / 3 + t;
152    T p3 = p / 3;
153    T q2 = q / 2;
154    T D  = p3 * p3 * p3 + q2 * q2;
155
156    if (D == 0 && p3 == 0)
157    {
158        x[0] = -r / 3;
159        x[1] = -r / 3;
160        x[2] = -r / 3;
161        return 1;
162    }
163
164    std::complex<T> u = std::pow (-q / 2 + std::sqrt (std::complex<T> (D)),
165                                  T (1) / T (3));
166
167    std::complex<T> v = -p / (T (3) * u);
168
169    const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
170                                                            // for long double
171    std::complex<T> y0 (u + v);
172
173    std::complex<T> y1 (-(u + v) / T (2) +
174                         (u - v) / T (2) * std::complex<T> (0, sqrt3));
175
176    std::complex<T> y2 (-(u + v) / T (2) -
177                         (u - v) / T (2) * std::complex<T> (0, sqrt3));
178
179    if (D > 0)
180    {
181        x[0] = y0.real() - r / 3;
182        return 1;
183    }
184    else if (D == 0)
185    {
186        x[0] = y0.real() - r / 3;
187        x[1] = y1.real() - r / 3;
188        return 2;
189    }
190    else
191    {
192        x[0] = y0.real() - r / 3;
193        x[1] = y1.real() - r / 3;
194        x[2] = y2.real() - r / 3;
195        return 3;
196    }
197}
198
199
200template <class T>
201int
202solveCubic (T a, T b, T c, T d, T x[3])
203{
204    if (a == 0)
205    {
206        return solveQuadratic (b, c, d, x);
207    }
208    else
209    {
210        return solveNormalizedCubic (b / a, c / a, d / a, x);
211    }
212}
213
214
215} // namespace Imath
216
217#endif
Note: See TracBrowser for help on using the repository browser.