OpenCores
URL https://opencores.org/ocsvn/altor32/altor32/trunk

Subversion Repositories altor32

[/] [altor32/] [trunk/] [gcc-x64/] [or1knd-elf/] [or1knd-elf/] [include/] [c++/] [4.8.0/] [tr1/] [ell_integral.tcc] - Blame information for rev 35

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 35 ultra_embe 
// Special functions -*- C++ -*-
2 
 
3 
// Copyright (C) 2006, 2007, 2008, 2009, 2010
4 
// Free Software Foundation, Inc.
5 
//
6 
// This file is part of the GNU ISO C++ Library.  This library is free
7 
// software; you can redistribute it and/or modify it under the
8 
// terms of the GNU General Public License as published by the
9 
// Free Software Foundation; either version 3, or (at your option)
10 
// any later version.
11 
//
12 
// This library is distributed in the hope that it will be useful,
13 
// but WITHOUT ANY WARRANTY; without even the implied warranty of
14 
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 
// GNU General Public License for more details.
16 
//
17 
// Under Section 7 of GPL version 3, you are granted additional
18 
// permissions described in the GCC Runtime Library Exception, version
19 
// 3.1, as published by the Free Software Foundation.
20 
 
21 
// You should have received a copy of the GNU General Public License and
22 
// a copy of the GCC Runtime Library Exception along with this program;
23 
// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
24 
// .
25 
 
26 
/** @file tr1/ell_integral.tcc
27 
 *  This is an internal header file, included by other library headers.
28 
 *  Do not attempt to use it directly. @headername{tr1/cmath}
29 
 */
30 
 
31 
//
32 
// ISO C++ 14882 TR1: 5.2  Special functions
33 
//
34 
 
35 
// Written by Edward Smith-Rowland based on:
36 
//   (1)  B. C. Carlson Numer. Math. 33, 1 (1979)
37 
//   (2)  B. C. Carlson, Special Functions of Applied Mathematics (1977)
38 
//   (3)  The Gnu Scientific Library, http://www.gnu.org/software/gsl
39 
//   (4)  Numerical Recipes in C, 2nd ed, by W. H. Press, S. A. Teukolsky,
40 
//        W. T. Vetterling, B. P. Flannery, Cambridge University Press
41 
//        (1992), pp. 261-269
42 
 
43 
#ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC
44 
#define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1
45 
 
46 
namespace std _GLIBCXX_VISIBILITY(default)
47 
{
48 
namespace tr1
49 
{
50 
  // [5.2] Special functions
51 
 
52 
  // Implementation-space details.
53 
  namespace __detail
54 
  {
55 
  _GLIBCXX_BEGIN_NAMESPACE_VERSION
56 
 
57 
    /**
58 
     *   @brief Return the Carlson elliptic function @f$ R_F(x,y,z) @f$
59 
     *          of the first kind.
60 
     *
61 
     *   The Carlson elliptic function of the first kind is defined by:
62 
     *   @f[
63 
     *       R_F(x,y,z) = \frac{1}{2} \int_0^\infty
64 
     *                 \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}}
65 
     *   @f]
66 
     *
67 
     *   @param  __x  The first of three symmetric arguments.
68 
     *   @param  __y  The second of three symmetric arguments.
69 
     *   @param  __z  The third of three symmetric arguments.
70 
     *   @return  The Carlson elliptic function of the first kind.
71 
     */
72 
    template
73 
    _Tp
74 
    __ellint_rf(const _Tp __x, const _Tp __y, const _Tp __z)
75 
    {
76 
      const _Tp __min = std::numeric_limits<_Tp>::min();
77 
      const _Tp __max = std::numeric_limits<_Tp>::max();
78 
      const _Tp __lolim = _Tp(5) * __min;
79 
      const _Tp __uplim = __max / _Tp(5);
80 
 
81 
      if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
82 
        std::__throw_domain_error(__N("Argument less than zero "
83 
                                      "in __ellint_rf."));
84 
      else if (__x + __y < __lolim || __x + __z < __lolim
85 
            || __y + __z < __lolim)
86 
        std::__throw_domain_error(__N("Argument too small in __ellint_rf"));
87 
      else
88 
        {
89 
          const _Tp __c0 = _Tp(1) / _Tp(4);
90 
          const _Tp __c1 = _Tp(1) / _Tp(24);
91 
          const _Tp __c2 = _Tp(1) / _Tp(10);
92 
          const _Tp __c3 = _Tp(3) / _Tp(44);
93 
          const _Tp __c4 = _Tp(1) / _Tp(14);
94 
 
95 
          _Tp __xn = __x;
96 
          _Tp __yn = __y;
97 
          _Tp __zn = __z;
98 
 
99 
          const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
100 
          const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6));
101 
          _Tp __mu;
102 
          _Tp __xndev, __yndev, __zndev;
103 
 
104 
          const unsigned int __max_iter = 100;
105 
          for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
106 
            {
107 
              __mu = (__xn + __yn + __zn) / _Tp(3);
108 
              __xndev = 2 - (__mu + __xn) / __mu;
109 
              __yndev = 2 - (__mu + __yn) / __mu;
110 
              __zndev = 2 - (__mu + __zn) / __mu;
111 
              _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
112 
              __epsilon = std::max(__epsilon, std::abs(__zndev));
113 
              if (__epsilon < __errtol)
114 
                break;
115 
              const _Tp __xnroot = std::sqrt(__xn);
116 
              const _Tp __ynroot = std::sqrt(__yn);
117 
              const _Tp __znroot = std::sqrt(__zn);
118 
              const _Tp __lambda = __xnroot * (__ynroot + __znroot)
119 
                                 + __ynroot * __znroot;
120 
              __xn = __c0 * (__xn + __lambda);
121 
              __yn = __c0 * (__yn + __lambda);
122 
              __zn = __c0 * (__zn + __lambda);
123 
            }
124 
 
125 
          const _Tp __e2 = __xndev * __yndev - __zndev * __zndev;
126 
          const _Tp __e3 = __xndev * __yndev * __zndev;
127 
          const _Tp __s  = _Tp(1) + (__c1 * __e2 - __c2 - __c3 * __e3) * __e2
128 
                   + __c4 * __e3;
129 
 
130 
          return __s / std::sqrt(__mu);
131 
        }
132 
    }
133 
 
134 
 
135 
    /**
136 
     *   @brief Return the complete elliptic integral of the first kind
137 
     *          @f$ K(k) @f$ by series expansion.
138 
     *
139 
     *   The complete elliptic integral of the first kind is defined as
140 
     *   @f[
141 
     *     K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
142 
     *                              {\sqrt{1 - k^2sin^2\theta}}
143 
     *   @f]
144 
     *
145 
     *   This routine is not bad as long as |k| is somewhat smaller than 1
146 
     *   but is not is good as the Carlson elliptic integral formulation.
147 
     *
148 
     *   @param  __k  The argument of the complete elliptic function.
149 
     *   @return  The complete elliptic function of the first kind.
150 
     */
151 
    template
152 
    _Tp
153 
    __comp_ellint_1_series(const _Tp __k)
154 
    {
155 
 
156 
      const _Tp __kk = __k * __k;
157 
 
158 
      _Tp __term = __kk / _Tp(4);
159 
      _Tp __sum = _Tp(1) + __term;
160 
 
161 
      const unsigned int __max_iter = 1000;
162 
      for (unsigned int __i = 2; __i < __max_iter; ++__i)
163 
        {
164 
          __term *= (2 * __i - 1) * __kk / (2 * __i);
165 
          if (__term < std::numeric_limits<_Tp>::epsilon())
166 
            break;
167 
          __sum += __term;
168 
        }
169 
 
170 
      return __numeric_constants<_Tp>::__pi_2() * __sum;
171 
    }
172 
 
173 
 
174 
    /**
175 
     *   @brief  Return the complete elliptic integral of the first kind
176 
     *           @f$ K(k) @f$ using the Carlson formulation.
177 
     *
178 
     *   The complete elliptic integral of the first kind is defined as
179 
     *   @f[
180 
     *     K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
181 
     *                                           {\sqrt{1 - k^2 sin^2\theta}}
182 
     *   @f]
183 
     *   where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the
184 
     *   first kind.
185 
     *
186 
     *   @param  __k  The argument of the complete elliptic function.
187 
     *   @return  The complete elliptic function of the first kind.
188 
     */
189 
    template
190 
    _Tp
191 
    __comp_ellint_1(const _Tp __k)
192 
    {
193 
 
194 
      if (__isnan(__k))
195 
        return std::numeric_limits<_Tp>::quiet_NaN();
196 
      else if (std::abs(__k) >= _Tp(1))
197 
        return std::numeric_limits<_Tp>::quiet_NaN();
198 
      else
199 
        return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1));
200 
    }
201 
 
202 
 
203 
    /**
204 
     *   @brief  Return the incomplete elliptic integral of the first kind
205 
     *           @f$ F(k,\phi) @f$ using the Carlson formulation.
206 
     *
207 
     *   The incomplete elliptic integral of the first kind is defined as
208 
     *   @f[
209 
     *     F(k,\phi) = \int_0^{\phi}\frac{d\theta}
210 
     *                                   {\sqrt{1 - k^2 sin^2\theta}}
211 
     *   @f]
212 
     *
213 
     *   @param  __k  The argument of the elliptic function.
214 
     *   @param  __phi  The integral limit argument of the elliptic function.
215 
     *   @return  The elliptic function of the first kind.
216 
     */
217 
    template
218 
    _Tp
219 
    __ellint_1(const _Tp __k, const _Tp __phi)
220 
    {
221 
 
222 
      if (__isnan(__k) || __isnan(__phi))
223 
        return std::numeric_limits<_Tp>::quiet_NaN();
224 
      else if (std::abs(__k) > _Tp(1))
225 
        std::__throw_domain_error(__N("Bad argument in __ellint_1."));
226 
      else
227 
        {
228 
          //  Reduce phi to -pi/2 < phi < +pi/2.
229 
          const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
230 
                                   + _Tp(0.5L));
231 
          const _Tp __phi_red = __phi
232 
                              - __n * __numeric_constants<_Tp>::__pi();
233 
 
234 
          const _Tp __s = std::sin(__phi_red);
235 
          const _Tp __c = std::cos(__phi_red);
236 
 
237 
          const _Tp __F = __s
238 
                        * __ellint_rf(__c * __c,
239 
                                _Tp(1) - __k * __k * __s * __s, _Tp(1));
240 
 
241 
          if (__n == 0)
242 
            return __F;
243 
          else
244 
            return __F + _Tp(2) * __n * __comp_ellint_1(__k);
245 
        }
246 
    }
247 
 
248 
 
249 
    /**
250 
     *   @brief Return the complete elliptic integral of the second kind
251 
     *          @f$ E(k) @f$ by series expansion.
252 
     *
253 
     *   The complete elliptic integral of the second kind is defined as
254 
     *   @f[
255 
     *     E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
256 
     *   @f]
257 
     *
258 
     *   This routine is not bad as long as |k| is somewhat smaller than 1
259 
     *   but is not is good as the Carlson elliptic integral formulation.
260 
     *
261 
     *   @param  __k  The argument of the complete elliptic function.
262 
     *   @return  The complete elliptic function of the second kind.
263 
     */
264 
    template
265 
    _Tp
266 
    __comp_ellint_2_series(const _Tp __k)
267 
    {
268 
 
269 
      const _Tp __kk = __k * __k;
270 
 
271 
      _Tp __term = __kk;
272 
      _Tp __sum = __term;
273 
 
274 
      const unsigned int __max_iter = 1000;
275 
      for (unsigned int __i = 2; __i < __max_iter; ++__i)
276 
        {
277 
          const _Tp __i2m = 2 * __i - 1;
278 
          const _Tp __i2 = 2 * __i;
279 
          __term *= __i2m * __i2m * __kk / (__i2 * __i2);
280 
          if (__term < std::numeric_limits<_Tp>::epsilon())
281 
            break;
282 
          __sum += __term / __i2m;
283 
        }
284 
 
285 
      return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum);
286 
    }
287 
 
288 
 
289 
    /**
290 
     *   @brief  Return the Carlson elliptic function of the second kind
291 
     *           @f$ R_D(x,y,z) = R_J(x,y,z,z) @f$ where
292 
     *           @f$ R_J(x,y,z,p) @f$ is the Carlson elliptic function
293 
     *           of the third kind.
294 
     *
295 
     *   The Carlson elliptic function of the second kind is defined by:
296 
     *   @f[
297 
     *       R_D(x,y,z) = \frac{3}{2} \int_0^\infty
298 
     *                 \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{3/2}}
299 
     *   @f]
300 
     *
301 
     *   Based on Carlson's algorithms:
302 
     *   -  B. C. Carlson Numer. Math. 33, 1 (1979)
303 
     *   -  B. C. Carlson, Special Functions of Applied Mathematics (1977)
304 
     *   -  Numerical Recipes in C, 2nd ed, pp. 261-269,
305 
     *      by Press, Teukolsky, Vetterling, Flannery (1992)
306 
     *
307 
     *   @param  __x  The first of two symmetric arguments.
308 
     *   @param  __y  The second of two symmetric arguments.
309 
     *   @param  __z  The third argument.
310 
     *   @return  The Carlson elliptic function of the second kind.
311 
     */
312 
    template
313 
    _Tp
314 
    __ellint_rd(const _Tp __x, const _Tp __y, const _Tp __z)
315 
    {
316 
      const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
317 
      const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
318 
      const _Tp __min = std::numeric_limits<_Tp>::min();
319 
      const _Tp __max = std::numeric_limits<_Tp>::max();
320 
      const _Tp __lolim = _Tp(2) / std::pow(__max, _Tp(2) / _Tp(3));
321 
      const _Tp __uplim = std::pow(_Tp(0.1L) * __errtol / __min, _Tp(2) / _Tp(3));
322 
 
323 
      if (__x < _Tp(0) || __y < _Tp(0))
324 
        std::__throw_domain_error(__N("Argument less than zero "
325 
                                      "in __ellint_rd."));
326 
      else if (__x + __y < __lolim || __z < __lolim)
327 
        std::__throw_domain_error(__N("Argument too small "
328 
                                      "in __ellint_rd."));
329 
      else
330 
        {
331 
          const _Tp __c0 = _Tp(1) / _Tp(4);
332 
          const _Tp __c1 = _Tp(3) / _Tp(14);
333 
          const _Tp __c2 = _Tp(1) / _Tp(6);
334 
          const _Tp __c3 = _Tp(9) / _Tp(22);
335 
          const _Tp __c4 = _Tp(3) / _Tp(26);
336 
 
337 
          _Tp __xn = __x;
338 
          _Tp __yn = __y;
339 
          _Tp __zn = __z;
340 
          _Tp __sigma = _Tp(0);
341 
          _Tp __power4 = _Tp(1);
342 
 
343 
          _Tp __mu;
344 
          _Tp __xndev, __yndev, __zndev;
345 
 
346 
          const unsigned int __max_iter = 100;
347 
          for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
348 
            {
349 
              __mu = (__xn + __yn + _Tp(3) * __zn) / _Tp(5);
350 
              __xndev = (__mu - __xn) / __mu;
351 
              __yndev = (__mu - __yn) / __mu;
352 
              __zndev = (__mu - __zn) / __mu;
353 
              _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
354 
              __epsilon = std::max(__epsilon, std::abs(__zndev));
355 
              if (__epsilon < __errtol)
356 
                break;
357 
              _Tp __xnroot = std::sqrt(__xn);
358 
              _Tp __ynroot = std::sqrt(__yn);
359 
              _Tp __znroot = std::sqrt(__zn);
360 
              _Tp __lambda = __xnroot * (__ynroot + __znroot)
361 
                           + __ynroot * __znroot;
362 
              __sigma += __power4 / (__znroot * (__zn + __lambda));
363 
              __power4 *= __c0;
364 
              __xn = __c0 * (__xn + __lambda);
365 
              __yn = __c0 * (__yn + __lambda);
366 
              __zn = __c0 * (__zn + __lambda);
367 
            }
368 
 
369 
          // Note: __ea is an SPU badname.
370 
          _Tp __eaa = __xndev * __yndev;
371 
          _Tp __eb = __zndev * __zndev;
372 
          _Tp __ec = __eaa - __eb;
373 
          _Tp __ed = __eaa - _Tp(6) * __eb;
374 
          _Tp __ef = __ed + __ec + __ec;
375 
          _Tp __s1 = __ed * (-__c1 + __c3 * __ed
376 
                                   / _Tp(3) - _Tp(3) * __c4 * __zndev * __ef
377 
                                   / _Tp(2));
378 
          _Tp __s2 = __zndev
379 
                   * (__c2 * __ef
380 
                    + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa));
381 
 
382 
          return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2)
383 
                                        / (__mu * std::sqrt(__mu));
384 
        }
385 
    }
386 
 
387 
 
388 
    /**
389 
     *   @brief  Return the complete elliptic integral of the second kind
390 
     *           @f$ E(k) @f$ using the Carlson formulation.
391 
     *
392 
     *   The complete elliptic integral of the second kind is defined as
393 
     *   @f[
394 
     *     E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
395 
     *   @f]
396 
     *
397 
     *   @param  __k  The argument of the complete elliptic function.
398 
     *   @return  The complete elliptic function of the second kind.
399 
     */
400 
    template
401 
    _Tp
402 
    __comp_ellint_2(const _Tp __k)
403 
    {
404 
 
405 
      if (__isnan(__k))
406 
        return std::numeric_limits<_Tp>::quiet_NaN();
407 
      else if (std::abs(__k) == 1)
408 
        return _Tp(1);
409 
      else if (std::abs(__k) > _Tp(1))
410 
        std::__throw_domain_error(__N("Bad argument in __comp_ellint_2."));
411 
      else
412 
        {
413 
          const _Tp __kk = __k * __k;
414 
 
415 
          return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
416 
               - __kk * __ellint_rd(_Tp(0), _Tp(1) - __kk, _Tp(1)) / _Tp(3);
417 
        }
418 
    }
419 
 
420 
 
421 
    /**
422 
     *   @brief  Return the incomplete elliptic integral of the second kind
423 
     *           @f$ E(k,\phi) @f$ using the Carlson formulation.
424 
     *
425 
     *   The incomplete elliptic integral of the second kind is defined as
426 
     *   @f[
427 
     *     E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta}
428 
     *   @f]
429 
     *
430 
     *   @param  __k  The argument of the elliptic function.
431 
     *   @param  __phi  The integral limit argument of the elliptic function.
432 
     *   @return  The elliptic function of the second kind.
433 
     */
434 
    template
435 
    _Tp
436 
    __ellint_2(const _Tp __k, const _Tp __phi)
437 
    {
438 
 
439 
      if (__isnan(__k) || __isnan(__phi))
440 
        return std::numeric_limits<_Tp>::quiet_NaN();
441 
      else if (std::abs(__k) > _Tp(1))
442 
        std::__throw_domain_error(__N("Bad argument in __ellint_2."));
443 
      else
444 
        {
445 
          //  Reduce phi to -pi/2 < phi < +pi/2.
446 
          const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
447 
                                   + _Tp(0.5L));
448 
          const _Tp __phi_red = __phi
449 
                              - __n * __numeric_constants<_Tp>::__pi();
450 
 
451 
          const _Tp __kk = __k * __k;
452 
          const _Tp __s = std::sin(__phi_red);
453 
          const _Tp __ss = __s * __s;
454 
          const _Tp __sss = __ss * __s;
455 
          const _Tp __c = std::cos(__phi_red);
456 
          const _Tp __cc = __c * __c;
457 
 
458 
          const _Tp __E = __s
459 
                        * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
460 
                        - __kk * __sss
461 
                        * __ellint_rd(__cc, _Tp(1) - __kk * __ss, _Tp(1))
462 
                        / _Tp(3);
463 
 
464 
          if (__n == 0)
465 
            return __E;
466 
          else
467 
            return __E + _Tp(2) * __n * __comp_ellint_2(__k);
468 
        }
469 
    }
470 
 
471 
 
472 
    /**
473 
     *   @brief  Return the Carlson elliptic function
474 
     *           @f$ R_C(x,y) = R_F(x,y,y) @f$ where @f$ R_F(x,y,z) @f$
475 
     *           is the Carlson elliptic function of the first kind.
476 
     *
477 
     *   The Carlson elliptic function is defined by:
478 
     *   @f[
479 
     *       R_C(x,y) = \frac{1}{2} \int_0^\infty
480 
     *                 \frac{dt}{(t + x)^{1/2}(t + y)}
481 
     *   @f]
482 
     *
483 
     *   Based on Carlson's algorithms:
484 
     *   -  B. C. Carlson Numer. Math. 33, 1 (1979)
485 
     *   -  B. C. Carlson, Special Functions of Applied Mathematics (1977)
486 
     *   -  Numerical Recipes in C, 2nd ed, pp. 261-269,
487 
     *      by Press, Teukolsky, Vetterling, Flannery (1992)
488 
     *
489 
     *   @param  __x  The first argument.
490 
     *   @param  __y  The second argument.
491 
     *   @return  The Carlson elliptic function.
492 
     */
493 
    template
494 
    _Tp
495 
    __ellint_rc(const _Tp __x, const _Tp __y)
496 
    {
497 
      const _Tp __min = std::numeric_limits<_Tp>::min();
498 
      const _Tp __max = std::numeric_limits<_Tp>::max();
499 
      const _Tp __lolim = _Tp(5) * __min;
500 
      const _Tp __uplim = __max / _Tp(5);
501 
 
502 
      if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim)
503 
        std::__throw_domain_error(__N("Argument less than zero "
504 
                                      "in __ellint_rc."));
505 
      else
506 
        {
507 
          const _Tp __c0 = _Tp(1) / _Tp(4);
508 
          const _Tp __c1 = _Tp(1) / _Tp(7);
509 
          const _Tp __c2 = _Tp(9) / _Tp(22);
510 
          const _Tp __c3 = _Tp(3) / _Tp(10);
511 
          const _Tp __c4 = _Tp(3) / _Tp(8);
512 
 
513 
          _Tp __xn = __x;
514 
          _Tp __yn = __y;
515 
 
516 
          const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
517 
          const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6));
518 
          _Tp __mu;
519 
          _Tp __sn;
520 
 
521 
          const unsigned int __max_iter = 100;
522 
          for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
523 
            {
524 
              __mu = (__xn + _Tp(2) * __yn) / _Tp(3);
525 
              __sn = (__yn + __mu) / __mu - _Tp(2);
526 
              if (std::abs(__sn) < __errtol)
527 
                break;
528 
              const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn)
529 
                             + __yn;
530 
              __xn = __c0 * (__xn + __lambda);
531 
              __yn = __c0 * (__yn + __lambda);
532 
            }
533 
 
534 
          _Tp __s = __sn * __sn
535 
                  * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2)));
536 
 
537 
          return (_Tp(1) + __s) / std::sqrt(__mu);
538 
        }
539 
    }
540 
 
541 
 
542 
    /**
543 
     *   @brief  Return the Carlson elliptic function @f$ R_J(x,y,z,p) @f$
544 
     *           of the third kind.
545 
     *
546 
     *   The Carlson elliptic function of the third kind is defined by:
547 
     *   @f[
548 
     *       R_J(x,y,z,p) = \frac{3}{2} \int_0^\infty
549 
     *       \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}(t + p)}
550 
     *   @f]
551 
     *
552 
     *   Based on Carlson's algorithms:
553 
     *   -  B. C. Carlson Numer. Math. 33, 1 (1979)
554 
     *   -  B. C. Carlson, Special Functions of Applied Mathematics (1977)
555 
     *   -  Numerical Recipes in C, 2nd ed, pp. 261-269,
556 
     *      by Press, Teukolsky, Vetterling, Flannery (1992)
557 
     *
558 
     *   @param  __x  The first of three symmetric arguments.
559 
     *   @param  __y  The second of three symmetric arguments.
560 
     *   @param  __z  The third of three symmetric arguments.
561 
     *   @param  __p  The fourth argument.
562 
     *   @return  The Carlson elliptic function of the fourth kind.
563 
     */
564 
    template
565 
    _Tp
566 
    __ellint_rj(const _Tp __x, const _Tp __y, const _Tp __z, const _Tp __p)
567 
    {
568 
      const _Tp __min = std::numeric_limits<_Tp>::min();
569 
      const _Tp __max = std::numeric_limits<_Tp>::max();
570 
      const _Tp __lolim = std::pow(_Tp(5) * __min, _Tp(1)/_Tp(3));
571 
      const _Tp __uplim = _Tp(0.3L)
572 
                        * std::pow(_Tp(0.2L) * __max, _Tp(1)/_Tp(3));
573 
 
574 
      if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0))
575 
        std::__throw_domain_error(__N("Argument less than zero "
576 
                                      "in __ellint_rj."));
577 
      else if (__x + __y < __lolim || __x + __z < __lolim
578 
            || __y + __z < __lolim || __p < __lolim)
579 
        std::__throw_domain_error(__N("Argument too small "
580 
                                      "in __ellint_rj"));
581 
      else
582 
        {
583 
          const _Tp __c0 = _Tp(1) / _Tp(4);
584 
          const _Tp __c1 = _Tp(3) / _Tp(14);
585 
          const _Tp __c2 = _Tp(1) / _Tp(3);
586 
          const _Tp __c3 = _Tp(3) / _Tp(22);
587 
          const _Tp __c4 = _Tp(3) / _Tp(26);
588 
 
589 
          _Tp __xn = __x;
590 
          _Tp __yn = __y;
591 
          _Tp __zn = __z;
592 
          _Tp __pn = __p;
593 
          _Tp __sigma = _Tp(0);
594 
          _Tp __power4 = _Tp(1);
595 
 
596 
          const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
597 
          const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
598 
 
599 
          _Tp __lambda, __mu;
600 
          _Tp __xndev, __yndev, __zndev, __pndev;
601 
 
602 
          const unsigned int __max_iter = 100;
603 
          for (unsigned int __iter = 0; __iter < __max_iter; ++__iter)
604 
            {
605 
              __mu = (__xn + __yn + __zn + _Tp(2) * __pn) / _Tp(5);
606 
              __xndev = (__mu - __xn) / __mu;
607 
              __yndev = (__mu - __yn) / __mu;
608 
              __zndev = (__mu - __zn) / __mu;
609 
              __pndev = (__mu - __pn) / __mu;
610 
              _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev));
611 
              __epsilon = std::max(__epsilon, std::abs(__zndev));
612 
              __epsilon = std::max(__epsilon, std::abs(__pndev));
613 
              if (__epsilon < __errtol)
614 
                break;
615 
              const _Tp __xnroot = std::sqrt(__xn);
616 
              const _Tp __ynroot = std::sqrt(__yn);
617 
              const _Tp __znroot = std::sqrt(__zn);
618 
              const _Tp __lambda = __xnroot * (__ynroot + __znroot)
619 
                                 + __ynroot * __znroot;
620 
              const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot)
621 
                                + __xnroot * __ynroot * __znroot;
622 
              const _Tp __alpha2 = __alpha1 * __alpha1;
623 
              const _Tp __beta = __pn * (__pn + __lambda)
624 
                                      * (__pn + __lambda);
625 
              __sigma += __power4 * __ellint_rc(__alpha2, __beta);
626 
              __power4 *= __c0;
627 
              __xn = __c0 * (__xn + __lambda);
628 
              __yn = __c0 * (__yn + __lambda);
629 
              __zn = __c0 * (__zn + __lambda);
630 
              __pn = __c0 * (__pn + __lambda);
631 
            }
632 
 
633 
          // Note: __ea is an SPU badname.
634 
          _Tp __eaa = __xndev * (__yndev + __zndev) + __yndev * __zndev;
635 
          _Tp __eb = __xndev * __yndev * __zndev;
636 
          _Tp __ec = __pndev * __pndev;
637 
          _Tp __e2 = __eaa - _Tp(3) * __ec;
638 
          _Tp __e3 = __eb + _Tp(2) * __pndev * (__eaa - __ec);
639 
          _Tp __s1 = _Tp(1) + __e2 * (-__c1 + _Tp(3) * __c3 * __e2 / _Tp(4)
640 
                            - _Tp(3) * __c4 * __e3 / _Tp(2));
641 
          _Tp __s2 = __eb * (__c2 / _Tp(2)
642 
                   + __pndev * (-__c3 - __c3 + __pndev * __c4));
643 
          _Tp __s3 = __pndev * __eaa * (__c2 - __pndev * __c3)
644 
                   - __c2 * __pndev * __ec;
645 
 
646 
          return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3)
647 
                                             / (__mu * std::sqrt(__mu));
648 
        }
649 
    }
650 
 
651 
 
652 
    /**
653 
     *   @brief Return the complete elliptic integral of the third kind
654 
     *          @f$ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) @f$ using the
655 
     *          Carlson formulation.
656 
     *
657 
     *   The complete elliptic integral of the third kind is defined as
658 
     *   @f[
659 
     *     \Pi(k,\nu) = \int_0^{\pi/2}
660 
     *                   \frac{d\theta}
661 
     *                 {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}}
662 
     *   @f]
663 
     *
664 
     *   @param  __k  The argument of the elliptic function.
665 
     *   @param  __nu  The second argument of the elliptic function.
666 
     *   @return  The complete elliptic function of the third kind.
667 
     */
668 
    template
669 
    _Tp
670 
    __comp_ellint_3(const _Tp __k, const _Tp __nu)
671 
    {
672 
 
673 
      if (__isnan(__k) || __isnan(__nu))
674 
        return std::numeric_limits<_Tp>::quiet_NaN();
675 
      else if (__nu == _Tp(1))
676 
        return std::numeric_limits<_Tp>::infinity();
677 
      else if (std::abs(__k) > _Tp(1))
678 
        std::__throw_domain_error(__N("Bad argument in __comp_ellint_3."));
679 
      else
680 
        {
681 
          const _Tp __kk = __k * __k;
682 
 
683 
          return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
684 
               - __nu
685 
               * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) + __nu)
686 
               / _Tp(3);
687 
        }
688 
    }
689 
 
690 
 
691 
    /**
692 
     *   @brief Return the incomplete elliptic integral of the third kind
693 
     *          @f$ \Pi(k,\nu,\phi) @f$ using the Carlson formulation.
694 
     *
695 
     *   The incomplete elliptic integral of the third kind is defined as
696 
     *   @f[
697 
     *     \Pi(k,\nu,\phi) = \int_0^{\phi}
698 
     *                       \frac{d\theta}
699 
     *                            {(1 - \nu \sin^2\theta)
700 
     *                             \sqrt{1 - k^2 \sin^2\theta}}
701 
     *   @f]
702 
     *
703 
     *   @param  __k  The argument of the elliptic function.
704 
     *   @param  __nu  The second argument of the elliptic function.
705 
     *   @param  __phi  The integral limit argument of the elliptic function.
706 
     *   @return  The elliptic function of the third kind.
707 
     */
708 
    template
709 
    _Tp
710 
    __ellint_3(const _Tp __k, const _Tp __nu, const _Tp __phi)
711 
    {
712 
 
713 
      if (__isnan(__k) || __isnan(__nu) || __isnan(__phi))
714 
        return std::numeric_limits<_Tp>::quiet_NaN();
715 
      else if (std::abs(__k) > _Tp(1))
716 
        std::__throw_domain_error(__N("Bad argument in __ellint_3."));
717 
      else
718 
        {
719 
          //  Reduce phi to -pi/2 < phi < +pi/2.
720 
          const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
721 
                                   + _Tp(0.5L));
722 
          const _Tp __phi_red = __phi
723 
                              - __n * __numeric_constants<_Tp>::__pi();
724 
 
725 
          const _Tp __kk = __k * __k;
726 
          const _Tp __s = std::sin(__phi_red);
727 
          const _Tp __ss = __s * __s;
728 
          const _Tp __sss = __ss * __s;
729 
          const _Tp __c = std::cos(__phi_red);
730 
          const _Tp __cc = __c * __c;
731 
 
732 
          const _Tp __Pi = __s
733 
                         * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
734 
                         - __nu * __sss
735 
                         * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1),
736 
                                       _Tp(1) + __nu * __ss) / _Tp(3);
737 
 
738 
          if (__n == 0)
739 
            return __Pi;
740 
          else
741 
            return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu);
742 
        }
743 
    }
744 
 
745 
  _GLIBCXX_END_NAMESPACE_VERSION
746 
  } // namespace std::tr1::__detail
747 
}
748 
}
749 
 
750 
#endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC
751 
 

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.