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

https://bitbucket.org/cabalistic/ogredeps/ · C++ Header · 218 lines · 115 code · 29 blank · 74 comment · 20 complexity · dbf465237ef05434d1699024844ecf3b MD5 · raw file

  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. #ifndef INCLUDED_IMATHROOTS_H
  35. #define INCLUDED_IMATHROOTS_H
  36. //---------------------------------------------------------------------
  37. //
  38. // Functions to solve linear, quadratic or cubic equations
  39. //
  40. //---------------------------------------------------------------------
  41. #include <ImathMath.h>
  42. #include <complex>
  43. namespace Imath {
  44. //--------------------------------------------------------------------------
  45. // Find the real solutions of a linear, quadratic or cubic equation:
  46. //
  47. // function equation solved
  48. //
  49. // solveLinear (a, b, x) a * x + b == 0
  50. // solveQuadratic (a, b, c, x) a * x*x + b * x + c == 0
  51. // solveNormalizedCubic (r, s, t, x) x*x*x + r * x*x + s * x + t == 0
  52. // solveCubic (a, b, c, d, x) a * x*x*x + b * x*x + c * x + d == 0
  53. //
  54. // Return value:
  55. //
  56. // 3 three real solutions, stored in x[0], x[1] and x[2]
  57. // 2 two real solutions, stored in x[0] and x[1]
  58. // 1 one real solution, stored in x[1]
  59. // 0 no real solutions
  60. // -1 all real numbers are solutions
  61. //
  62. // Notes:
  63. //
  64. // * It is possible that an equation has real solutions, but that the
  65. // solutions (or some intermediate result) are not representable.
  66. // In this case, either some of the solutions returned are invalid
  67. // (nan or infinity), or, if floating-point exceptions have been
  68. // enabled with Iex::mathExcOn(), an Iex::MathExc exception is
  69. // thrown.
  70. //
  71. // * Cubic equations are solved using Cardano's Formula; even though
  72. // only real solutions are produced, some intermediate results are
  73. // complex (std::complex<T>).
  74. //
  75. //--------------------------------------------------------------------------
  76. template <class T> int solveLinear (T a, T b, T &x);
  77. template <class T> int solveQuadratic (T a, T b, T c, T x[2]);
  78. template <class T> int solveNormalizedCubic (T r, T s, T t, T x[3]);
  79. template <class T> int solveCubic (T a, T b, T c, T d, T x[3]);
  80. //---------------
  81. // Implementation
  82. //---------------
  83. template <class T>
  84. int
  85. solveLinear (T a, T b, T &x)
  86. {
  87. if (a != 0)
  88. {
  89. x = -b / a;
  90. return 1;
  91. }
  92. else if (b != 0)
  93. {
  94. return 0;
  95. }
  96. else
  97. {
  98. return -1;
  99. }
  100. }
  101. template <class T>
  102. int
  103. solveQuadratic (T a, T b, T c, T x[2])
  104. {
  105. if (a == 0)
  106. {
  107. return solveLinear (b, c, x[0]);
  108. }
  109. else
  110. {
  111. T D = b * b - 4 * a * c;
  112. if (D > 0)
  113. {
  114. T s = Math<T>::sqrt (D);
  115. x[0] = (-b + s) / (2 * a);
  116. x[1] = (-b - s) / (2 * a);
  117. return 2;
  118. }
  119. if (D == 0)
  120. {
  121. x[0] = -b / (2 * a);
  122. return 1;
  123. }
  124. else
  125. {
  126. return 0;
  127. }
  128. }
  129. }
  130. template <class T>
  131. int
  132. solveNormalizedCubic (T r, T s, T t, T x[3])
  133. {
  134. T p = (3 * s - r * r) / 3;
  135. T q = 2 * r * r * r / 27 - r * s / 3 + t;
  136. T p3 = p / 3;
  137. T q2 = q / 2;
  138. T D = p3 * p3 * p3 + q2 * q2;
  139. if (D == 0 && p3 == 0)
  140. {
  141. x[0] = -r / 3;
  142. x[1] = -r / 3;
  143. x[2] = -r / 3;
  144. return 1;
  145. }
  146. std::complex<T> u = std::pow (-q / 2 + std::sqrt (std::complex<T> (D)),
  147. T (1) / T (3));
  148. std::complex<T> v = -p / (T (3) * u);
  149. const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
  150. // for long double
  151. std::complex<T> y0 (u + v);
  152. std::complex<T> y1 (-(u + v) / T (2) +
  153. (u - v) / T (2) * std::complex<T> (0, sqrt3));
  154. std::complex<T> y2 (-(u + v) / T (2) -
  155. (u - v) / T (2) * std::complex<T> (0, sqrt3));
  156. if (D > 0)
  157. {
  158. x[0] = y0.real() - r / 3;
  159. return 1;
  160. }
  161. else if (D == 0)
  162. {
  163. x[0] = y0.real() - r / 3;
  164. x[1] = y1.real() - r / 3;
  165. return 2;
  166. }
  167. else
  168. {
  169. x[0] = y0.real() - r / 3;
  170. x[1] = y1.real() - r / 3;
  171. x[2] = y2.real() - r / 3;
  172. return 3;
  173. }
  174. }
  175. template <class T>
  176. int
  177. solveCubic (T a, T b, T c, T d, T x[3])
  178. {
  179. if (a == 0)
  180. {
  181. return solveQuadratic (b, c, d, x);
  182. }
  183. else
  184. {
  185. return solveNormalizedCubic (b / a, c / a, d / a, x);
  186. }
  187. }
  188. } // namespace Imath
  189. #endif