/src/FreeImage/Source/OpenEXR/Imath/ImathRoots.h

  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
130	    x[0] = (-b + s) / (2 * a);
131	    x[1] = (-b - s) / (2 * a);
132	    return 2;
133	}
134	if (D == 0)
135	{
136	    x[0] = -b / (2 * a);
137	    return 1;
138	}
139	else
140	{
141	    return 0;
142	}
143    }
144}
145
146
147template <class T>
148int
149solveNormalizedCubic (T r, T s, T t, T x[3])
150{
151    T p  = (3 * s - r * r) / 3;
152    T q  = 2 * r * r * r / 27 - r * s / 3 + t;
153    T p3 = p / 3;
154    T q2 = q / 2;
155    T D  = p3 * p3 * p3 + q2 * q2;
156
157    if (D == 0 && p3 == 0)
158    {
159	x[0] = -r / 3;
160	x[1] = -r / 3;
161	x[2] = -r / 3;
162	return 1;
163    }
164
165    std::complex<T> u = std::pow (-q / 2 + std::sqrt (std::complex<T> (D)),
166				  T (1) / T (3));
167
168    std::complex<T> v = -p / (T (3) * u);
169
170    const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
171							    // for long double
172    std::complex<T> y0 (u + v);
173
174    std::complex<T> y1 (-(u + v) / T (2) +
175			 (u - v) / T (2) * std::complex<T> (0, sqrt3));
176
177    std::complex<T> y2 (-(u + v) / T (2) -
178			 (u - v) / T (2) * std::complex<T> (0, sqrt3));
179
180    if (D > 0)
181    {
182	x[0] = y0.real() - r / 3;
183	return 1;
184    }
185    else if (D == 0)
186    {
187	x[0] = y0.real() - r / 3;
188	x[1] = y1.real() - r / 3;
189	return 2;
190    }
191    else
192    {
193	x[0] = y0.real() - r / 3;
194	x[1] = y1.real() - r / 3;
195	x[2] = y2.real() - r / 3;
196	return 3;
197    }
198}
199
200
201template <class T>
202int
203solveCubic (T a, T b, T c, T d, T x[3])
204{
205    if (a == 0)
206    {
207	return solveQuadratic (b, c, d, x);
208    }
209    else
210    {
211	return solveNormalizedCubic (b / a, c / a, d / a, x);
212    }
213}
214
215
216} // namespace Imath
217
218#endif