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