| 1 |
// Special functions -*- C++ -*- |
| 2 |
|
| 3 |
// Copyright (C) 2006-2021 Free Software Foundation, Inc. |
| 4 |
// |
| 5 |
// This file is part of the GNU ISO C++ Library. This library is free |
| 6 |
// software; you can redistribute it and/or modify it under the |
| 7 |
// terms of the GNU General Public License as published by the |
| 8 |
// Free Software Foundation; either version 3, or (at your option) |
| 9 |
// any later version. |
| 10 |
// |
| 11 |
// This library is distributed in the hope that it will be useful, |
| 12 |
// but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 13 |
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 14 |
// GNU General Public License for more details. |
| 15 |
// |
| 16 |
// Under Section 7 of GPL version 3, you are granted additional |
| 17 |
// permissions described in the GCC Runtime Library Exception, version |
| 18 |
// 3.1, as published by the Free Software Foundation. |
| 19 |
|
| 20 |
// You should have received a copy of the GNU General Public License and |
| 21 |
// a copy of the GCC Runtime Library Exception along with this program; |
| 22 |
// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
| 23 |
// <http://www.gnu.org/licenses/>. |
| 24 |
|
| 25 |
/** @file tr1/legendre_function.tcc |
| 26 |
* This is an internal header file, included by other library headers. |
| 27 |
* Do not attempt to use it directly. @headername{tr1/cmath} |
| 28 |
*/ |
| 29 |
|
| 30 |
// |
| 31 |
// ISO C++ 14882 TR1: 5.2 Special functions |
| 32 |
// |
| 33 |
|
| 34 |
// Written by Edward Smith-Rowland based on: |
| 35 |
// (1) Handbook of Mathematical Functions, |
| 36 |
// ed. Milton Abramowitz and Irene A. Stegun, |
| 37 |
// Dover Publications, |
| 38 |
// Section 8, pp. 331-341 |
| 39 |
// (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl |
| 40 |
// (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, |
| 41 |
// W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), |
| 42 |
// 2nd ed, pp. 252-254 |
| 43 |
|
| 44 |
#ifndef _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC |
| 45 |
#define _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC 1 |
| 46 |
|
| 47 |
#include <tr1/special_function_util.h> |
| 48 |
|
| 49 |
namespace std _GLIBCXX_VISIBILITY(default) |
| 50 |
{ |
| 51 |
_GLIBCXX_BEGIN_NAMESPACE_VERSION |
| 52 |
|
| 53 |
#if _GLIBCXX_USE_STD_SPEC_FUNCS |
| 54 |
# define _GLIBCXX_MATH_NS ::std |
| 55 |
#elif defined(_GLIBCXX_TR1_CMATH) |
| 56 |
namespace tr1 |
| 57 |
{ |
| 58 |
# define _GLIBCXX_MATH_NS ::std::tr1 |
| 59 |
#else |
| 60 |
# error do not include this header directly, use <cmath> or <tr1/cmath> |
| 61 |
#endif |
| 62 |
// [5.2] Special functions |
| 63 |
|
| 64 |
// Implementation-space details. |
| 65 |
namespace __detail |
| 66 |
{ |
| 67 |
/** |
| 68 |
* @brief Return the Legendre polynomial by recursion on degree |
| 69 |
* @f$ l @f$. |
| 70 |
* |
| 71 |
* The Legendre function of @f$ l @f$ and @f$ x @f$, |
| 72 |
* @f$ P_l(x) @f$, is defined by: |
| 73 |
* @f[ |
| 74 |
* P_l(x) = \frac{1}{2^l l!}\frac{d^l}{dx^l}(x^2 - 1)^{l} |
| 75 |
* @f] |
| 76 |
* |
| 77 |
* @param l The degree of the Legendre polynomial. @f$l >= 0@f$. |
| 78 |
* @param x The argument of the Legendre polynomial. @f$|x| <= 1@f$. |
| 79 |
*/ |
| 80 |
template<typename _Tp> |
| 81 |
_Tp |
| 82 |
__poly_legendre_p(unsigned int __l, _Tp __x) |
| 83 |
{ |
| 84 |
|
| 85 |
if (__isnan(__x)) |
| 86 |
return std::numeric_limits<_Tp>::quiet_NaN(); |
| 87 |
else if (__x == +_Tp(1)) |
| 88 |
return +_Tp(1); |
| 89 |
else if (__x == -_Tp(1)) |
| 90 |
return (__l % 2 == 1 ? -_Tp(1) : +_Tp(1)); |
| 91 |
else |
| 92 |
{ |
| 93 |
_Tp __p_lm2 = _Tp(1); |
| 94 |
if (__l == 0) |
| 95 |
return __p_lm2; |
| 96 |
|
| 97 |
_Tp __p_lm1 = __x; |
| 98 |
if (__l == 1) |
| 99 |
return __p_lm1; |
| 100 |
|
| 101 |
_Tp __p_l = 0; |
| 102 |
for (unsigned int __ll = 2; __ll <= __l; ++__ll) |
| 103 |
{ |
| 104 |
// This arrangement is supposed to be better for roundoff |
| 105 |
// protection, Arfken, 2nd Ed, Eq 12.17a. |
| 106 |
__p_l = _Tp(2) * __x * __p_lm1 - __p_lm2 |
| 107 |
- (__x * __p_lm1 - __p_lm2) / _Tp(__ll); |
| 108 |
__p_lm2 = __p_lm1; |
| 109 |
__p_lm1 = __p_l; |
| 110 |
} |
| 111 |
|
| 112 |
return __p_l; |
| 113 |
} |
| 114 |
} |
| 115 |
|
| 116 |
|
| 117 |
/** |
| 118 |
* @brief Return the associated Legendre function by recursion |
| 119 |
* on @f$ l @f$. |
| 120 |
* |
| 121 |
* The associated Legendre function is derived from the Legendre function |
| 122 |
* @f$ P_l(x) @f$ by the Rodrigues formula: |
| 123 |
* @f[ |
| 124 |
* P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x) |
| 125 |
* @f] |
| 126 |
* @note @f$ P_l^m(x) = 0 @f$ if @f$ m > l @f$. |
| 127 |
* |
| 128 |
* @param l The degree of the associated Legendre function. |
| 129 |
* @f$ l >= 0 @f$. |
| 130 |
* @param m The order of the associated Legendre function. |
| 131 |
* @param x The argument of the associated Legendre function. |
| 132 |
* @f$ |x| <= 1 @f$. |
| 133 |
* @param phase The phase of the associated Legendre function. |
| 134 |
* Use -1 for the Condon-Shortley phase convention. |
| 135 |
*/ |
| 136 |
template<typename _Tp> |
| 137 |
_Tp |
| 138 |
__assoc_legendre_p(unsigned int __l, unsigned int __m, _Tp __x, |
| 139 |
_Tp __phase = _Tp(+1)) |
| 140 |
{ |
| 141 |
|
| 142 |
if (__m > __l) |
| 143 |
return _Tp(0); |
| 144 |
else if (__isnan(__x)) |
| 145 |
return std::numeric_limits<_Tp>::quiet_NaN(); |
| 146 |
else if (__m == 0) |
| 147 |
return __poly_legendre_p(__l, __x); |
| 148 |
else |
| 149 |
{ |
| 150 |
_Tp __p_mm = _Tp(1); |
| 151 |
if (__m > 0) |
| 152 |
{ |
| 153 |
// Two square roots seem more accurate more of the time |
| 154 |
// than just one. |
| 155 |
_Tp __root = std::sqrt(_Tp(1) - __x) * std::sqrt(_Tp(1) + __x); |
| 156 |
_Tp __fact = _Tp(1); |
| 157 |
for (unsigned int __i = 1; __i <= __m; ++__i) |
| 158 |
{ |
| 159 |
__p_mm *= __phase * __fact * __root; |
| 160 |
__fact += _Tp(2); |
| 161 |
} |
| 162 |
} |
| 163 |
if (__l == __m) |
| 164 |
return __p_mm; |
| 165 |
|
| 166 |
_Tp __p_mp1m = _Tp(2 * __m + 1) * __x * __p_mm; |
| 167 |
if (__l == __m + 1) |
| 168 |
return __p_mp1m; |
| 169 |
|
| 170 |
_Tp __p_lm2m = __p_mm; |
| 171 |
_Tp __P_lm1m = __p_mp1m; |
| 172 |
_Tp __p_lm = _Tp(0); |
| 173 |
for (unsigned int __j = __m + 2; __j <= __l; ++__j) |
| 174 |
{ |
| 175 |
__p_lm = (_Tp(2 * __j - 1) * __x * __P_lm1m |
| 176 |
- _Tp(__j + __m - 1) * __p_lm2m) / _Tp(__j - __m); |
| 177 |
__p_lm2m = __P_lm1m; |
| 178 |
__P_lm1m = __p_lm; |
| 179 |
} |
| 180 |
|
| 181 |
return __p_lm; |
| 182 |
} |
| 183 |
} |
| 184 |
|
| 185 |
|
| 186 |
/** |
| 187 |
* @brief Return the spherical associated Legendre function. |
| 188 |
* |
| 189 |
* The spherical associated Legendre function of @f$ l @f$, @f$ m @f$, |
| 190 |
* and @f$ \theta @f$ is defined as @f$ Y_l^m(\theta,0) @f$ where |
| 191 |
* @f[ |
| 192 |
* Y_l^m(\theta,\phi) = (-1)^m[\frac{(2l+1)}{4\pi} |
| 193 |
* \frac{(l-m)!}{(l+m)!}] |
| 194 |
* P_l^m(\cos\theta) \exp^{im\phi} |
| 195 |
* @f] |
| 196 |
* is the spherical harmonic function and @f$ P_l^m(x) @f$ is the |
| 197 |
* associated Legendre function. |
| 198 |
* |
| 199 |
* This function differs from the associated Legendre function by |
| 200 |
* argument (@f$x = \cos(\theta)@f$) and by a normalization factor |
| 201 |
* but this factor is rather large for large @f$ l @f$ and @f$ m @f$ |
| 202 |
* and so this function is stable for larger differences of @f$ l @f$ |
| 203 |
* and @f$ m @f$. |
| 204 |
* @note Unlike the case for __assoc_legendre_p the Condon-Shortley |
| 205 |
* phase factor @f$ (-1)^m @f$ is present here. |
| 206 |
* @note @f$ Y_l^m(\theta) = 0 @f$ if @f$ m > l @f$. |
| 207 |
* |
| 208 |
* @param l The degree of the spherical associated Legendre function. |
| 209 |
* @f$ l >= 0 @f$. |
| 210 |
* @param m The order of the spherical associated Legendre function. |
| 211 |
* @param theta The radian angle argument of the spherical associated |
| 212 |
* Legendre function. |
| 213 |
*/ |
| 214 |
template <typename _Tp> |
| 215 |
_Tp |
| 216 |
__sph_legendre(unsigned int __l, unsigned int __m, _Tp __theta) |
| 217 |
{ |
| 218 |
if (__isnan(__theta)) |
| 219 |
return std::numeric_limits<_Tp>::quiet_NaN(); |
| 220 |
|
| 221 |
const _Tp __x = std::cos(__theta); |
| 222 |
|
| 223 |
if (__m > __l) |
| 224 |
return _Tp(0); |
| 225 |
else if (__m == 0) |
| 226 |
{ |
| 227 |
_Tp __P = __poly_legendre_p(__l, __x); |
| 228 |
_Tp __fact = std::sqrt(_Tp(2 * __l + 1) |
| 229 |
/ (_Tp(4) * __numeric_constants<_Tp>::__pi())); |
| 230 |
__P *= __fact; |
| 231 |
return __P; |
| 232 |
} |
| 233 |
else if (__x == _Tp(1) || __x == -_Tp(1)) |
| 234 |
{ |
| 235 |
// m > 0 here |
| 236 |
return _Tp(0); |
| 237 |
} |
| 238 |
else |
| 239 |
{ |
| 240 |
// m > 0 and |x| < 1 here |
| 241 |
|
| 242 |
// Starting value for recursion. |
| 243 |
// Y_m^m(x) = sqrt( (2m+1)/(4pi m) gamma(m+1/2)/gamma(m) ) |
| 244 |
// (-1)^m (1-x^2)^(m/2) / pi^(1/4) |
| 245 |
const _Tp __sgn = ( __m % 2 == 1 ? -_Tp(1) : _Tp(1)); |
| 246 |
const _Tp __y_mp1m_factor = __x * std::sqrt(_Tp(2 * __m + 3)); |
| 247 |
#if _GLIBCXX_USE_C99_MATH_TR1 |
| 248 |
const _Tp __lncirc = _GLIBCXX_MATH_NS::log1p(-__x * __x); |
| 249 |
#else |
| 250 |
const _Tp __lncirc = std::log(_Tp(1) - __x * __x); |
| 251 |
#endif |
| 252 |
// Gamma(m+1/2) / Gamma(m) |
| 253 |
#if _GLIBCXX_USE_C99_MATH_TR1 |
| 254 |
const _Tp __lnpoch = _GLIBCXX_MATH_NS::lgamma(_Tp(__m + _Tp(0.5L))) |
| 255 |
- _GLIBCXX_MATH_NS::lgamma(_Tp(__m)); |
| 256 |
#else |
| 257 |
const _Tp __lnpoch = __log_gamma(_Tp(__m + _Tp(0.5L))) |
| 258 |
- __log_gamma(_Tp(__m)); |
| 259 |
#endif |
| 260 |
const _Tp __lnpre_val = |
| 261 |
-_Tp(0.25L) * __numeric_constants<_Tp>::__lnpi() |
| 262 |
+ _Tp(0.5L) * (__lnpoch + __m * __lncirc); |
| 263 |
const _Tp __sr = std::sqrt((_Tp(2) + _Tp(1) / __m) |
| 264 |
/ (_Tp(4) * __numeric_constants<_Tp>::__pi())); |
| 265 |
_Tp __y_mm = __sgn * __sr * std::exp(__lnpre_val); |
| 266 |
_Tp __y_mp1m = __y_mp1m_factor * __y_mm; |
| 267 |
|
| 268 |
if (__l == __m) |
| 269 |
return __y_mm; |
| 270 |
else if (__l == __m + 1) |
| 271 |
return __y_mp1m; |
| 272 |
else |
| 273 |
{ |
| 274 |
_Tp __y_lm = _Tp(0); |
| 275 |
|
| 276 |
// Compute Y_l^m, l > m+1, upward recursion on l. |
| 277 |
for (unsigned int __ll = __m + 2; __ll <= __l; ++__ll) |
| 278 |
{ |
| 279 |
const _Tp __rat1 = _Tp(__ll - __m) / _Tp(__ll + __m); |
| 280 |
const _Tp __rat2 = _Tp(__ll - __m - 1) / _Tp(__ll + __m - 1); |
| 281 |
const _Tp __fact1 = std::sqrt(__rat1 * _Tp(2 * __ll + 1) |
| 282 |
* _Tp(2 * __ll - 1)); |
| 283 |
const _Tp __fact2 = std::sqrt(__rat1 * __rat2 * _Tp(2 * __ll + 1) |
| 284 |
/ _Tp(2 * __ll - 3)); |
| 285 |
__y_lm = (__x * __y_mp1m * __fact1 |
| 286 |
- (__ll + __m - 1) * __y_mm * __fact2) / _Tp(__ll - __m); |
| 287 |
__y_mm = __y_mp1m; |
| 288 |
__y_mp1m = __y_lm; |
| 289 |
} |
| 290 |
|
| 291 |
return __y_lm; |
| 292 |
} |
| 293 |
} |
| 294 |
} |
| 295 |
} // namespace __detail |
| 296 |
#undef _GLIBCXX_MATH_NS |
| 297 |
#if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) |
| 298 |
} // namespace tr1 |
| 299 |
#endif |
| 300 |
|
| 301 |
_GLIBCXX_END_NAMESPACE_VERSION |
| 302 |
} |
| 303 |
|
| 304 |
#endif // _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC |
