1// Special functions -*- C++ -*- 2 3// Copyright (C) 2006-2015 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 "special_function_util.h" 48 49namespace std _GLIBCXX_VISIBILITY(default) 50{ 51namespace tr1 52{ 53 // [5.2] Special functions 54 55 // Implementation-space details. 56 namespace __detail 57 { 58 _GLIBCXX_BEGIN_NAMESPACE_VERSION 59 60 /** 61 * @brief Return the Legendre polynomial by recursion on order 62 * @f$ l @f$. 63 * 64 * The Legendre function of @f$ l @f$ and @f$ x @f$, 65 * @f$ P_l(x) @f$, is defined by: 66 * @f[ 67 * P_l(x) = \frac{1}{2^l l!}\frac{d^l}{dx^l}(x^2 - 1)^{l} 68 * @f] 69 * 70 * @param l The order of the Legendre polynomial. @f$l >= 0@f$. 71 * @param x The argument of the Legendre polynomial. @f$|x| <= 1@f$. 72 */ 73 template<typename _Tp> 74 _Tp 75 __poly_legendre_p(unsigned int __l, _Tp __x) 76 { 77 78 if ((__x < _Tp(-1)) || (__x > _Tp(+1))) 79 std::__throw_domain_error(__N("Argument out of range" 80 " in __poly_legendre_p.")); 81 else if (__isnan(__x)) 82 return std::numeric_limits<_Tp>::quiet_NaN(); 83 else if (__x == +_Tp(1)) 84 return +_Tp(1); 85 else if (__x == -_Tp(1)) 86 return (__l % 2 == 1 ? -_Tp(1) : +_Tp(1)); 87 else 88 { 89 _Tp __p_lm2 = _Tp(1); 90 if (__l == 0) 91 return __p_lm2; 92 93 _Tp __p_lm1 = __x; 94 if (__l == 1) 95 return __p_lm1; 96 97 _Tp __p_l = 0; 98 for (unsigned int __ll = 2; __ll <= __l; ++__ll) 99 { 100 // This arrangement is supposed to be better for roundoff 101 // protection, Arfken, 2nd Ed, Eq 12.17a. 102 __p_l = _Tp(2) * __x * __p_lm1 - __p_lm2 103 - (__x * __p_lm1 - __p_lm2) / _Tp(__ll); 104 __p_lm2 = __p_lm1; 105 __p_lm1 = __p_l; 106 } 107 108 return __p_l; 109 } 110 } 111 112 113 /** 114 * @brief Return the associated Legendre function by recursion 115 * on @f$ l @f$. 116 * 117 * The associated Legendre function is derived from the Legendre function 118 * @f$ P_l(x) @f$ by the Rodrigues formula: 119 * @f[ 120 * P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x) 121 * @f] 122 * 123 * @param l The order of the associated Legendre function. 124 * @f$ l >= 0 @f$. 125 * @param m The order of the associated Legendre function. 126 * @f$ m <= l @f$. 127 * @param x The argument of the associated Legendre function. 128 * @f$ |x| <= 1 @f$. 129 */ 130 template<typename _Tp> 131 _Tp 132 __assoc_legendre_p(unsigned int __l, unsigned int __m, _Tp __x) 133 { 134 135 if (__x < _Tp(-1) || __x > _Tp(+1)) 136 std::__throw_domain_error(__N("Argument out of range" 137 " in __assoc_legendre_p.")); 138 else if (__m > __l) 139 std::__throw_domain_error(__N("Degree out of range" 140 " in __assoc_legendre_p.")); 141 else if (__isnan(__x)) 142 return std::numeric_limits<_Tp>::quiet_NaN(); 143 else if (__m == 0) 144 return __poly_legendre_p(__l, __x); 145 else 146 { 147 _Tp __p_mm = _Tp(1); 148 if (__m > 0) 149 { 150 // Two square roots seem more accurate more of the time 151 // than just one. 152 _Tp __root = std::sqrt(_Tp(1) - __x) * std::sqrt(_Tp(1) + __x); 153 _Tp __fact = _Tp(1); 154 for (unsigned int __i = 1; __i <= __m; ++__i) 155 { 156 __p_mm *= -__fact * __root; 157 __fact += _Tp(2); 158 } 159 } 160 if (__l == __m) 161 return __p_mm; 162 163 _Tp __p_mp1m = _Tp(2 * __m + 1) * __x * __p_mm; 164 if (__l == __m + 1) 165 return __p_mp1m; 166 167 _Tp __p_lm2m = __p_mm; 168 _Tp __P_lm1m = __p_mp1m; 169 _Tp __p_lm = _Tp(0); 170 for (unsigned int __j = __m + 2; __j <= __l; ++__j) 171 { 172 __p_lm = (_Tp(2 * __j - 1) * __x * __P_lm1m 173 - _Tp(__j + __m - 1) * __p_lm2m) / _Tp(__j - __m); 174 __p_lm2m = __P_lm1m; 175 __P_lm1m = __p_lm; 176 } 177 178 return __p_lm; 179 } 180 } 181 182 183 /** 184 * @brief Return the spherical associated Legendre function. 185 * 186 * The spherical associated Legendre function of @f$ l @f$, @f$ m @f$, 187 * and @f$ \theta @f$ is defined as @f$ Y_l^m(\theta,0) @f$ where 188 * @f[ 189 * Y_l^m(\theta,\phi) = (-1)^m[\frac{(2l+1)}{4\pi} 190 * \frac{(l-m)!}{(l+m)!}] 191 * P_l^m(\cos\theta) \exp^{im\phi} 192 * @f] 193 * is the spherical harmonic function and @f$ P_l^m(x) @f$ is the 194 * associated Legendre function. 195 * 196 * This function differs from the associated Legendre function by 197 * argument (@f$x = \cos(\theta)@f$) and by a normalization factor 198 * but this factor is rather large for large @f$ l @f$ and @f$ m @f$ 199 * and so this function is stable for larger differences of @f$ l @f$ 200 * and @f$ m @f$. 201 * 202 * @param l The order of the spherical associated Legendre function. 203 * @f$ l >= 0 @f$. 204 * @param m The order of the spherical associated Legendre function. 205 * @f$ m <= l @f$. 206 * @param theta The radian angle argument of the spherical associated 207 * Legendre function. 208 */ 209 template <typename _Tp> 210 _Tp 211 __sph_legendre(unsigned int __l, unsigned int __m, _Tp __theta) 212 { 213 if (__isnan(__theta)) 214 return std::numeric_limits<_Tp>::quiet_NaN(); 215 216 const _Tp __x = std::cos(__theta); 217 218 if (__l < __m) 219 { 220 std::__throw_domain_error(__N("Bad argument " 221 "in __sph_legendre.")); 222 } 223 else if (__m == 0) 224 { 225 _Tp __P = __poly_legendre_p(__l, __x); 226 _Tp __fact = std::sqrt(_Tp(2 * __l + 1) 227 / (_Tp(4) * __numeric_constants<_Tp>::__pi())); 228 __P *= __fact; 229 return __P; 230 } 231 else if (__x == _Tp(1) || __x == -_Tp(1)) 232 { 233 // m > 0 here 234 return _Tp(0); 235 } 236 else 237 { 238 // m > 0 and |x| < 1 here 239 240 // Starting value for recursion. 241 // Y_m^m(x) = sqrt( (2m+1)/(4pi m) gamma(m+1/2)/gamma(m) ) 242 // (-1)^m (1-x^2)^(m/2) / pi^(1/4) 243 const _Tp __sgn = ( __m % 2 == 1 ? -_Tp(1) : _Tp(1)); 244 const _Tp __y_mp1m_factor = __x * std::sqrt(_Tp(2 * __m + 3)); 245#if _GLIBCXX_USE_C99_MATH_TR1 246 const _Tp __lncirc = std::tr1::log1p(-__x * __x); 247#else 248 const _Tp __lncirc = std::log(_Tp(1) - __x * __x); 249#endif 250 // Gamma(m+1/2) / Gamma(m) 251#if _GLIBCXX_USE_C99_MATH_TR1 252 const _Tp __lnpoch = std::tr1::lgamma(_Tp(__m + _Tp(0.5L))) 253 - std::tr1::lgamma(_Tp(__m)); 254#else 255 const _Tp __lnpoch = __log_gamma(_Tp(__m + _Tp(0.5L))) 256 - __log_gamma(_Tp(__m)); 257#endif 258 const _Tp __lnpre_val = 259 -_Tp(0.25L) * __numeric_constants<_Tp>::__lnpi() 260 + _Tp(0.5L) * (__lnpoch + __m * __lncirc); 261 _Tp __sr = std::sqrt((_Tp(2) + _Tp(1) / __m) 262 / (_Tp(4) * __numeric_constants<_Tp>::__pi())); 263 _Tp __y_mm = __sgn * __sr * std::exp(__lnpre_val); 264 _Tp __y_mp1m = __y_mp1m_factor * __y_mm; 265 266 if (__l == __m) 267 { 268 return __y_mm; 269 } 270 else if (__l == __m + 1) 271 { 272 return __y_mp1m; 273 } 274 else 275 { 276 _Tp __y_lm = _Tp(0); 277 278 // Compute Y_l^m, l > m+1, upward recursion on l. 279 for ( int __ll = __m + 2; __ll <= __l; ++__ll) 280 { 281 const _Tp __rat1 = _Tp(__ll - __m) / _Tp(__ll + __m); 282 const _Tp __rat2 = _Tp(__ll - __m - 1) / _Tp(__ll + __m - 1); 283 const _Tp __fact1 = std::sqrt(__rat1 * _Tp(2 * __ll + 1) 284 * _Tp(2 * __ll - 1)); 285 const _Tp __fact2 = std::sqrt(__rat1 * __rat2 * _Tp(2 * __ll + 1) 286 / _Tp(2 * __ll - 3)); 287 __y_lm = (__x * __y_mp1m * __fact1 288 - (__ll + __m - 1) * __y_mm * __fact2) / _Tp(__ll - __m); 289 __y_mm = __y_mp1m; 290 __y_mp1m = __y_lm; 291 } 292 293 return __y_lm; 294 } 295 } 296 } 297 298 _GLIBCXX_END_NAMESPACE_VERSION 299 } // namespace std::tr1::__detail 300} 301} 302 303#endif // _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC 304