← Back to branch summary

~ubuntu-branches/ubuntu/saucy/openexr/saucy

« back to all changes in this revision

Viewing changes to Imath/ImathRoots.h

  • Committer: Bazaar Package Importer
  • Author(s): Adeodato Simó
  • Date: 2008-03-24 23:00:21 UTC
  • mfrom: (3.1.2 lenny)
  • Revision ID: james.westby@ubuntu.com-20080324230021-gnofz9mnvcj1xlv3
Tags: 1.6.1-3
Disable (hopefully temporarily) the test suite on arm and ia64.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Imath/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