← Back to branch summary

~cosme/ubuntu/precise/freeimage/freeimage-3.15.1

« back to all changes in this revision

Viewing changes to Source/OpenEXR/Imath/ImathRoots.h

  • Committer: Stefano Rivera
  • Date: 2010-07-24 15:35:51 UTC
  • mto: This revision was merged to the branch mainline in revision 5.
  • Revision ID: stefanor@ubuntu.com-20100724153551-6s3fth1653huk31a
Tags: upstream-3.13.1
Import upstream version 3.31.1

Show diffs side-by-side

added added

removed removed

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