← Back to branch summary

~ubuntu-branches/ubuntu/saucy/merkaartor/saucy

« back to all changes in this revision

Viewing changes to include/builtin-boost/boost/math/special_functions/ellint_3.hpp

Tags: upstream-0.15.3+svn20934
Import upstream version 0.15.3+svn20934

Show diffs side-by-side

added added

removed removed

Lines of Context:
include/builtin-boost/boost/math/special_functions/ellint_3.hpp
 
1
//  Copyright (c) 2006 Xiaogang Zhang
 
2
//  Copyright (c) 2006 John Maddock
 
3
//  Use, modification and distribution are subject to the
 
4
//  Boost Software License, Version 1.0. (See accompanying file
 
5
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 
6
//
 
7
//  History:
 
8
//  XZ wrote the original of this file as part of the Google
 
9
//  Summer of Code 2006.  JM modified it to fit into the
 
10
//  Boost.Math conceptual framework better, and to correctly
 
11
//  handle the various corner cases.
 
12
//
 
13
 
 
14
#ifndef BOOST_MATH_ELLINT_3_HPP
 
15
#define BOOST_MATH_ELLINT_3_HPP
 
16
 
 
17
#ifdef _MSC_VER
 
18
#pragma once
 
19
#endif
 
20
 
 
21
#include <boost/math/special_functions/ellint_rf.hpp>
 
22
#include <boost/math/special_functions/ellint_rj.hpp>
 
23
#include <boost/math/special_functions/ellint_1.hpp>
 
24
#include <boost/math/special_functions/ellint_2.hpp>
 
25
#include <boost/math/special_functions/log1p.hpp>
 
26
#include <boost/math/constants/constants.hpp>
 
27
#include <boost/math/policies/error_handling.hpp>
 
28
#include <boost/math/tools/workaround.hpp>
 
29
 
 
30
// Elliptic integrals (complete and incomplete) of the third kind
 
31
// Carlson, Numerische Mathematik, vol 33, 1 (1979)
 
32
 
 
33
namespace boost { namespace math { 
 
34
   
 
35
namespace detail{
 
36
 
 
37
template <typename T, typename Policy>
 
38
T ellint_pi_imp(T v, T k, T vc, const Policy& pol);
 
39
 
 
40
// Elliptic integral (Legendre form) of the third kind
 
41
template <typename T, typename Policy>
 
42
T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
 
43
{
 
44
    // Note vc = 1-v presumably without cancellation error.
 
45
    T value, x, y, z, p, t;
 
46
 
 
47
    BOOST_MATH_STD_USING
 
48
    using namespace boost::math::tools;
 
49
    using namespace boost::math::constants;
 
50
 
 
51
    static const char* function = "boost::math::ellint_3<%1%>(%1%,%1%,%1%)";
 
52
 
 
53
    if (abs(k) > 1)
 
54
    {
 
55
       return policies::raise_domain_error<T>(function,
 
56
            "Got k = %1%, function requires |k| <= 1", k, pol);
 
57
    }
 
58
 
 
59
    T sphi = sin(fabs(phi));
 
60
 
 
61
    if(v > 1 / (sphi * sphi))
 
62
    {
 
63
        // Complex result is a domain error:
 
64
       return policies::raise_domain_error<T>(function,
 
65
            "Got v = %1%, but result is complex for v > 1 / sin^2(phi)", v, pol);
 
66
    }
 
67
 
 
68
    // Special cases first:
 
69
    if(v == 0)
 
70
    {
 
71
       // A&S 17.7.18 & 19
 
72
       return (k == 0) ? phi : ellint_f_imp(phi, k, pol);
 
73
    }
 
74
    if(phi == constants::pi<T>() / 2)
 
75
    {
 
76
       // Have to filter this case out before the next
 
77
       // special case, otherwise we might get an infinity from
 
78
       // tan(phi).
 
79
       // Also note that since we can't represent PI/2 exactly
 
80
       // in a T, this is a bit of a guess as to the users true
 
81
       // intent...
 
82
       //
 
83
       return ellint_pi_imp(v, k, vc, pol);
 
84
    }
 
85
    if(k == 0)
 
86
    {
 
87
       // A&S 17.7.20:
 
88
       if(v < 1)
 
89
       {
 
90
          T vcr = sqrt(vc);
 
91
          return atan(vcr * tan(phi)) / vcr;
 
92
       }
 
93
       else if(v == 1)
 
94
       {
 
95
          return tan(phi);
 
96
       }
 
97
       else
 
98
       {
 
99
          // v > 1:
 
100
          T vcr = sqrt(-vc);
 
101
          T arg = vcr * tan(phi);
 
102
          return (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * vcr);
 
103
       }
 
104
    }
 
105
 
 
106
    if(v < 0)
 
107
    {
 
108
       //
 
109
       // If we don't shift to 0 <= v <= 1 we get
 
110
       // cancellation errors later on.  Use
 
111
       // A&S 17.7.15/16 to shift to v > 0:
 
112
       //
 
113
       T k2 = k * k;
 
114
       T N = (k2 - v) / (1 - v);
 
115
       T Nm1 = (1 - k2) / (1 - v);
 
116
       T p2 = sqrt(-v * (k2 - v) / (1 - v));
 
117
       T delta = sqrt(1 - k2 * sphi * sphi);
 
118
       T result = ellint_pi_imp(N, phi, k, Nm1, pol);
 
119
 
 
120
       result *= sqrt(Nm1 * (1 - k2 / N));
 
121
       result += ellint_f_imp(phi, k, pol) * k2 / p2;
 
122
       result += atan((p2/2) * sin(2 * phi) / delta);
 
123
       result /= sqrt((1 - v) * (1 - k2 / v));
 
124
       return result;
 
125
    }
 
126
#if 0  // disabled but retained for future reference: see below.
 
127
    if(v > 1)
 
128
    {
 
129
       //
 
130
       // If v > 1 we can use the identity in A&S 17.7.7/8
 
131
       // to shift to 0 <= v <= 1.  Unfortunately this
 
132
       // identity appears only to function correctly when
 
133
       // 0 <= phi <= pi/2, but it's when phi is outside that
 
134
       // range that we really need it: That's when
 
135
       // Carlson's formula fails, and what's more the periodicity
 
136
       // reduction used below on phi doesn't work when v > 1.
 
137
       //
 
138
       // So we're stuck... the code is archived here in case
 
139
       // some bright spart can figure out the fix.
 
140
       //
 
141
       T k2 = k * k;
 
142
       T N = k2 / v;
 
143
       T Nm1 = (v - k2) / v;
 
144
       T p1 = sqrt((-vc) * (1 - k2 / v));
 
145
       T delta = sqrt(1 - k2 * sphi * sphi);
 
146
       //
 
147
       // These next two terms have a large amount of cancellation
 
148
       // so it's not clear if this relation is useable even if
 
149
       // the issues with phi > pi/2 can be fixed:
 
150
       //
 
151
       T result = -ellint_pi_imp(N, phi, k, Nm1);
 
152
       result += ellint_f_imp(phi, k);
 
153
       //
 
154
       // This log term gives the complex result when
 
155
       //     n > 1/sin^2(phi)
 
156
       // However that case is dealt with as an error above, 
 
157
       // so we should always get a real result here:
 
158
       //
 
159
       result += log((delta + p1 * tan(phi)) / (delta - p1 * tan(phi))) / (2 * p1);
 
160
       return result;
 
161
    }
 
162
#endif
 
163
 
 
164
    // Carlson's algorithm works only for |phi| <= pi/2,
 
165
    // use the integrand's periodicity to normalize phi
 
166
    //
 
167
    // Xiaogang's original code used a cast to long long here
 
168
    // but that fails if T has more digits than a long long,
 
169
    // so rewritten to use fmod instead:
 
170
    //
 
171
    if(fabs(phi) > 1 / tools::epsilon<T>())
 
172
    {
 
173
       if(v > 1)
 
174
          return policies::raise_domain_error<T>(
 
175
            function,
 
176
            "Got v = %1%, but this is only supported for 0 <= phi <= pi/2", v, pol);
 
177
       //  
 
178
       // Phi is so large that phi%pi is necessarily zero (or garbage),
 
179
       // just return the second part of the duplication formula:
 
180
       //
 
181
       value = 2 * fabs(phi) * ellint_pi_imp(v, k, vc, pol) / constants::pi<T>();
 
182
    }
 
183
    else
 
184
    {
 
185
       T rphi = boost::math::tools::fmod_workaround(fabs(phi), constants::pi<T>() / 2);
 
186
       T m = floor((2 * fabs(phi)) / constants::pi<T>());
 
187
       int sign = 1;
 
188
       if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
 
189
       {
 
190
          m += 1;
 
191
          sign = -1;
 
192
          rphi = constants::pi<T>() / 2 - rphi;
 
193
       }
 
194
 
 
195
       if((m > 0) && (v > 1))
 
196
       {
 
197
          //
 
198
          // The region with v > 1 and phi outside [0, pi/2] is
 
199
          // currently unsupported:
 
200
          //
 
201
          return policies::raise_domain_error<T>(
 
202
            function,
 
203
            "Got v = %1%, but this is only supported for 0 <= phi <= pi/2", v, pol);
 
204
       }  
 
205
       T sinp = sin(rphi);
 
206
       T cosp = cos(rphi);
 
207
       x = cosp * cosp;
 
208
       t = sinp * sinp;
 
209
       y = 1 - k * k * t;
 
210
       z = 1;
 
211
       if(v * t < 0.5)
 
212
           p = 1 - v * t;
 
213
       else
 
214
           p = x + vc * t;
 
215
       value = sign * sinp * (ellint_rf_imp(x, y, z, pol) + v * t * ellint_rj_imp(x, y, z, p, pol) / 3);
 
216
       if(m > 0)
 
217
         value += m * ellint_pi_imp(v, k, vc, pol);
 
218
    }
 
219
 
 
220
    if (phi < 0)
 
221
    {
 
222
        value = -value;    // odd function
 
223
    }
 
224
    return value;
 
225
}
 
226
 
 
227
// Complete elliptic integral (Legendre form) of the third kind
 
228
template <typename T, typename Policy>
 
229
T ellint_pi_imp(T v, T k, T vc, const Policy& pol)
 
230
{
 
231
    // Note arg vc = 1-v, possibly without cancellation errors
 
232
    BOOST_MATH_STD_USING
 
233
    using namespace boost::math::tools;
 
234
 
 
235
    static const char* function = "boost::math::ellint_pi<%1%>(%1%,%1%)";
 
236
 
 
237
    if (abs(k) >= 1)
 
238
    {
 
239
       return policies::raise_domain_error<T>(function,
 
240
            "Got k = %1%, function requires |k| <= 1", k, pol);
 
241
    }
 
242
    if(vc <= 0)
 
243
    {
 
244
       // Result is complex:
 
245
       return policies::raise_domain_error<T>(function,
 
246
            "Got v = %1%, function requires v < 1", v, pol);
 
247
    }
 
248
 
 
249
    if(v == 0)
 
250
    {
 
251
       return (k == 0) ? boost::math::constants::pi<T>() / 2 : ellint_k_imp(k, pol);
 
252
    }
 
253
 
 
254
    if(v < 0)
 
255
    {
 
256
       T k2 = k * k;
 
257
       T N = (k2 - v) / (1 - v);
 
258
       T Nm1 = (1 - k2) / (1 - v);
 
259
       T p2 = sqrt(-v * (k2 - v) / (1 - v));
 
260
 
 
261
       T result = boost::math::detail::ellint_pi_imp(N, k, Nm1, pol);
 
262
 
 
263
       result *= sqrt(Nm1 * (1 - k2 / N));
 
264
       result += ellint_k_imp(k, pol) * k2 / p2;
 
265
       result /= sqrt((1 - v) * (1 - k2 / v));
 
266
       return result;
 
267
    }
 
268
 
 
269
    T x = 0;
 
270
    T y = 1 - k * k;
 
271
    T z = 1;
 
272
    T p = vc;
 
273
    T value = ellint_rf_imp(x, y, z, pol) + v * ellint_rj_imp(x, y, z, p, pol) / 3;
 
274
 
 
275
    return value;
 
276
}
 
277
 
 
278
template <class T1, class T2, class T3>
 
279
inline typename tools::promote_args<T1, T2, T3>::type ellint_3(T1 k, T2 v, T3 phi, const mpl::false_&)
 
280
{
 
281
   return boost::math::ellint_3(k, v, phi, policies::policy<>());
 
282
}
 
283
 
 
284
template <class T1, class T2, class Policy>
 
285
inline typename tools::promote_args<T1, T2>::type ellint_3(T1 k, T2 v, const Policy& pol, const mpl::true_&)
 
286
{
 
287
   typedef typename tools::promote_args<T1, T2>::type result_type;
 
288
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
 
289
   return policies::checked_narrowing_cast<result_type, Policy>(
 
290
      detail::ellint_pi_imp(
 
291
         static_cast<value_type>(v), 
 
292
         static_cast<value_type>(k),
 
293
         static_cast<value_type>(1-v),
 
294
         pol), "boost::math::ellint_3<%1%>(%1%,%1%)");
 
295
}
 
296
 
 
297
} // namespace detail
 
298
 
 
299
template <class T1, class T2, class T3, class Policy>
 
300
inline typename tools::promote_args<T1, T2, T3>::type ellint_3(T1 k, T2 v, T3 phi, const Policy& pol)
 
301
{
 
302
   typedef typename tools::promote_args<T1, T2, T3>::type result_type;
 
303
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
 
304
   return policies::checked_narrowing_cast<result_type, Policy>(
 
305
      detail::ellint_pi_imp(
 
306
         static_cast<value_type>(v), 
 
307
         static_cast<value_type>(phi), 
 
308
         static_cast<value_type>(k),
 
309
         static_cast<value_type>(1-v),
 
310
         pol), "boost::math::ellint_3<%1%>(%1%,%1%,%1%)");
 
311
}
 
312
 
 
313
template <class T1, class T2, class T3>
 
314
typename detail::ellint_3_result<T1, T2, T3>::type ellint_3(T1 k, T2 v, T3 phi)
 
315
{
 
316
   typedef typename policies::is_policy<T3>::type tag_type;
 
317
   return detail::ellint_3(k, v, phi, tag_type());
 
318
}
 
319
 
 
320
template <class T1, class T2>
 
321
inline typename tools::promote_args<T1, T2>::type ellint_3(T1 k, T2 v)
 
322
{
 
323
   return ellint_3(k, v, policies::policy<>());
 
324
}
 
325
 
 
326
}} // namespaces
 
327
 
 
328
#endif // BOOST_MATH_ELLINT_3_HPP
 
329