PBRT
|
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