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