PBRT
/home/felix/UBC/projects/AdaptiveLightfieldSampling/pbrt_v2/src/3rdparty/ilmbase-1.0.2/ImathRoots.h
00001 
00002 //
00003 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
00004 // Digital Ltd. LLC
00005 // 
00006 // All rights reserved.
00007 // 
00008 // Redistribution and use in source and binary forms, with or without
00009 // modification, are permitted provided that the following conditions are
00010 // met:
00011 // *       Redistributions of source code must retain the above copyright
00012 // notice, this list of conditions and the following disclaimer.
00013 // *       Redistributions in binary form must reproduce the above
00014 // copyright notice, this list of conditions and the following disclaimer
00015 // in the documentation and/or other materials provided with the
00016 // distribution.
00017 // *       Neither the name of Industrial Light & Magic nor the names of
00018 // its contributors may be used to endorse or promote products derived
00019 // from this software without specific prior written permission. 
00020 // 
00021 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00022 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00023 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00024 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
00025 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
00027 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
00028 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
00029 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
00030 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
00031 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00032 //
00034 
00035 
00036 
00037 #ifndef INCLUDED_IMATHROOTS_H
00038 #define INCLUDED_IMATHROOTS_H
00039 
00040 //---------------------------------------------------------------------
00041 //
00042 //      Functions to solve linear, quadratic or cubic equations
00043 //
00044 //---------------------------------------------------------------------
00045 
00046 #include <ImathMath.h>
00047 #include <complex>
00048 
00049 namespace Imath {
00050 
00051 //--------------------------------------------------------------------------
00052 // Find the real solutions of a linear, quadratic or cubic equation:
00053 //
00054 //      function                                   equation solved
00055 //
00056 //   solveLinear (a, b, x)                                    a * x + b == 0
00057 //   solveQuadratic (a, b, c, x)                    a * x*x + b * x + c == 0
00058 //   solveNormalizedCubic (r, s, t, x)      x*x*x + r * x*x + s * x + t == 0
00059 //   solveCubic (a, b, c, d, x)         a * x*x*x + b * x*x + c * x + d == 0
00060 //
00061 // Return value:
00062 //
00063 //       3      three real solutions, stored in x[0], x[1] and x[2]
00064 //       2      two real solutions, stored in x[0] and x[1]
00065 //       1      one real solution, stored in x[1]
00066 //       0      no real solutions
00067 //      -1      all real numbers are solutions
00068 //
00069 // Notes:
00070 //
00071 //    * It is possible that an equation has real solutions, but that the
00072 //      solutions (or some intermediate result) are not representable.
00073 //      In this case, either some of the solutions returned are invalid
00074 //      (nan or infinity), or, if floating-point exceptions have been
00075 //      enabled with Iex::mathExcOn(), an Iex::MathExc exception is
00076 //      thrown.
00077 //
00078 //    * Cubic equations are solved using Cardano's Formula; even though
00079 //      only real solutions are produced, some intermediate results are
00080 //      complex (std::complex<T>).
00081 //
00082 //--------------------------------------------------------------------------
00083 
00084 template <class T> int  solveLinear (T a, T b, T &x);
00085 template <class T> int  solveQuadratic (T a, T b, T c, T x[2]);
00086 template <class T> int  solveNormalizedCubic (T r, T s, T t, T x[3]);
00087 template <class T> int  solveCubic (T a, T b, T c, T d, T x[3]);
00088 
00089 
00090 //---------------
00091 // Implementation
00092 //---------------
00093 
00094 template <class T>
00095 int
00096 solveLinear (T a, T b, T &x)
00097 {
00098     if (a != 0)
00099     {
00100         x = -b / a;
00101         return 1;
00102     }
00103     else if (b != 0)
00104     {
00105         return 0;
00106     }
00107     else
00108     {
00109         return -1;
00110     }
00111 }
00112 
00113 
00114 template <class T>
00115 int
00116 solveQuadratic (T a, T b, T c, T x[2])
00117 {
00118     if (a == 0)
00119     {
00120         return solveLinear (b, c, x[0]);
00121     }
00122     else
00123     {
00124         T D = b * b - 4 * a * c;
00125 
00126         if (D > 0)
00127         {
00128             T s = Math<T>::sqrt (D);
00129 
00130             x[0] = (-b + s) / (2 * a);
00131             x[1] = (-b - s) / (2 * a);
00132             return 2;
00133         }
00134         if (D == 0)
00135         {
00136             x[0] = -b / (2 * a);
00137             return 1;
00138         }
00139         else
00140         {
00141             return 0;
00142         }
00143     }
00144 }
00145 
00146 
00147 template <class T>
00148 int
00149 solveNormalizedCubic (T r, T s, T t, T x[3])
00150 {
00151     T p  = (3 * s - r * r) / 3;
00152     T q  = 2 * r * r * r / 27 - r * s / 3 + t;
00153     T p3 = p / 3;
00154     T q2 = q / 2;
00155     T D  = p3 * p3 * p3 + q2 * q2;
00156 
00157     if (D == 0 && p3 == 0)
00158     {
00159         x[0] = -r / 3;
00160         x[1] = -r / 3;
00161         x[2] = -r / 3;
00162         return 1;
00163     }
00164 
00165     std::complex<T> u = std::pow (-q / 2 + std::sqrt (std::complex<T> (D)),
00166                                   T (1) / T (3));
00167 
00168     std::complex<T> v = -p / (T (3) * u);
00169 
00170     const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
00171                                                             // for long double
00172     std::complex<T> y0 (u + v);
00173 
00174     std::complex<T> y1 (-(u + v) / T (2) +
00175                          (u - v) / T (2) * std::complex<T> (0, sqrt3));
00176 
00177     std::complex<T> y2 (-(u + v) / T (2) -
00178                          (u - v) / T (2) * std::complex<T> (0, sqrt3));
00179 
00180     if (D > 0)
00181     {
00182         x[0] = y0.real() - r / 3;
00183         return 1;
00184     }
00185     else if (D == 0)
00186     {
00187         x[0] = y0.real() - r / 3;
00188         x[1] = y1.real() - r / 3;
00189         return 2;
00190     }
00191     else
00192     {
00193         x[0] = y0.real() - r / 3;
00194         x[1] = y1.real() - r / 3;
00195         x[2] = y2.real() - r / 3;
00196         return 3;
00197     }
00198 }
00199 
00200 
00201 template <class T>
00202 int
00203 solveCubic (T a, T b, T c, T d, T x[3])
00204 {
00205     if (a == 0)
00206     {
00207         return solveQuadratic (b, c, d, x);
00208     }
00209     else
00210     {
00211         return solveNormalizedCubic (b / a, c / a, d / a, x);
00212     }
00213 }
00214 
00215 
00216 } // namespace Imath
00217 
00218 #endif