1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006-2024 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 305
Indexes created Tue Mar 03 05:31:39 UTC 2026