diff options
Diffstat (limited to 'libstdc++-v3/include/bits/sf_gegenbauer.tcc')
| -rw-r--r-- | libstdc++-v3/include/bits/sf_gegenbauer.tcc | 143 |
1 files changed, 125 insertions, 18 deletions
diff --git a/libstdc++-v3/include/bits/sf_gegenbauer.tcc b/libstdc++-v3/include/bits/sf_gegenbauer.tcc index ef3be01adec..4c12dc30c06 100644 --- a/libstdc++-v3/include/bits/sf_gegenbauer.tcc +++ b/libstdc++-v3/include/bits/sf_gegenbauer.tcc @@ -55,36 +55,143 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION * @tparam _Talpha The real type of the order * @tparam _Tp The real type of the argument * @param __n The non-negative integral degree - * @param __alpha The real order + * @param __alpha1 The real order * @param __x The real argument */ template<typename _Tp> _Tp - __gegenbauer_poly(unsigned int __n, _Tp __alpha, _Tp __x) + __gegenbauer_poly(unsigned int __n, _Tp __alpha1, _Tp __x) { const auto _S_NaN = __gnu_cxx::__quiet_NaN(__x); - if (__isnan(__alpha) || __isnan(__x)) + if (__isnan(__alpha1) || __isnan(__x)) return _S_NaN; - auto _C0 = _Tp{1}; + auto __C_nm2 = _Tp{1}; if (__n == 0) - return _C0; + return __C_nm2; - auto _C1 = _Tp{2} * __alpha * __x; + auto __C_nm1 = _Tp{2} * __alpha1 * __x; if (__n == 1) - return _C1; - - auto _Cn = _Tp{0}; - for (unsigned int __nn = 2; __nn <= __n; ++__nn) - { - _Cn = (_Tp{2} * (_Tp{__nn} - _Tp{1} + __alpha) * __x * _C1 - - (_Tp{__nn} - _Tp{2} + _Tp{2} * __alpha) * _C0) - / _Tp(__nn); - _C0 = _C1; - _C1 = _Cn; - } - return _Cn; + return __C_nm1; + + auto __C_n = (_Tp{2} * (_Tp{1} + __alpha1) * __x * __C_nm1 + - _Tp{2} * __alpha1 * __C_nm2) / _Tp(2); + for (unsigned int __nn = 3; __nn <= __n; ++__nn) + { + __C_nm2 = __C_nm1; + __C_nm1 = __C_n; + __C_n = (_Tp{2} * (_Tp(__nn) - _Tp{1} + __alpha1) * __x * __C_nm1 + - (_Tp(__nn) - _Tp{2} + _Tp{2} * __alpha1) * __C_nm2) + / _Tp(__nn); + } + return __C_n; + } + + /** + * Return a vector containing the zeros of the Gegenbauer or ultraspherical + * polynomial @f$ C_n^{(\alpha)}@f$. + */ + template<typename _Tp> + std::vector<__gnu_cxx::__quadrature_point_t<_Tp>> + __gegenbauer_zeros(unsigned int __n, _Tp __alpha1) + { + const auto _S_eps = __gnu_cxx::__epsilon(__alpha1); + const unsigned int _S_maxit = 1000u; + std::vector<__gnu_cxx::__quadrature_point_t<_Tp>> __pt(__n); + + _Tp __z; + _Tp __w; + for (auto __i = 1u; __i <= __n; ++__i) + { + if (__i == 1) + { + auto __an = __alpha1 / __n; + auto __an2 = __an * __an; + auto __r1 = (1.0 + __alpha1) * (2.78 / (4.0 + __n * __n) + + 0.768 * __an / __n); + auto __r2 = 1.0 + 1.48 * __an + 0.96 * __an + 1.282 * __an2; + __z = 1.0 - __r1 / __r2; + } + else if (__i == 2) + { + auto __r1 = (4.1 + __alpha1) + / ((1.0 + __alpha1) * (1.0 + 0.156 * __alpha1)); + auto __r2 = 1.0 + + 0.06 * (__n - 8.0) * (1.0 + 0.12 * __alpha1) / __n; + auto __r3 = 1.0 + + 0.012 * __alpha1 * (1.0 + 0.25 * std::abs(__alpha1)) / __n; + __z -= (1.0 - __z) * __r1 * __r2 * __r3; + } + else if (__i == 3) + { + auto __r1 = (1.67 + 0.28 * __alpha1) / (1.0 + 0.37 * __alpha1); + auto __r2 = 1.0 + 0.22 * (__n - 8.0) / __n; + auto __r3 = 1.0 + 8.0 *__alpha1 / ((6.28 + __alpha1) * __n * __n); + __z -= (__pt[0].__zero - __z) * __r1 * __r2 * __r3; + } + else if (__i == __n - 1) + { + auto __r1 = (1.0 + 0.235 * __alpha1) / (0.766 + 0.119 * __alpha1); + auto __r2 = 1.0 / (1.0 + 0.639 * (__n - 4.0) + / (1.0 + 0.71 * (__n - 4.0))); + auto __r3 = 1.0 / (1.0 + 20.0 * __alpha1 + / ((7.5 + __alpha1) * __n * __n)); + __z += (__z - __pt[__n - 4].__zero) * __r1 * __r2 * __r3; + } + else if (__i == __n) + { + auto __r1 = (1.0 + 0.37 * __alpha1) / (1.67 + 0.28 * __alpha1); + auto __r2 = 1.0 / (1.0 + 0.22 * (__n - 8.0) / __n); + auto __r3 = 1.0 / (1.0 + 8.0 * __alpha1 + / ((6.28 + __alpha1) * __n * __n)); + __z += (__z - __pt[__n - 3].__zero) * __r1 * __r2 * __r3; + } + else + __z = 3.0 * __pt[__i - 2].__zero + - 3.0 * __pt[__i - 3].__zero + __pt[__i - 4].__zero; + + auto __2alpha = _Tp{2} * __alpha1; + for (auto __its = 1u; __its <= _S_maxit; ++__its) + { + auto __temp = _Tp{2} + __2alpha; + auto __C1 = (__temp * __z) / _Tp{2}; + auto __C2 = _Tp{1}; + for (auto __j = 2u; __j <= __n; ++__j) + { + auto __C3 = __C2; + __C2 = __C1; + __temp = _Tp{2} * __j + __2alpha; + auto __a = _Tp{2} * __j * (__j + __2alpha) + * (__temp - _Tp{2}); + auto __b = (__temp - _Tp{1}) + * __temp * (__temp - _Tp{2}) * __z; + auto __c = _Tp{2} * (__j - 1 + __alpha1) + * (__j - 1 + __alpha1) * __temp; + __C1 = (__b * __C2 - __c * __C3) / __a; + } + auto __Cp = (__n * (-__temp * __z) * __C1 + + _Tp{2} * (__n + __alpha1) * (__n + __alpha1) * __C2) + / (__temp * (_Tp{1} - __z * __z)); + auto __z1 = __z; + __z = __z1 - __C1 / __Cp; + if (std::abs(__z - __z1) <= _S_eps) + { + __w = std::exp(std::lgamma(__alpha1 + _Tp(__n)) + + std::lgamma(__alpha1 + _Tp(__n)) + - std::lgamma(_Tp(__n + 1)) + - std::lgamma(_Tp(__n + 1) + __2alpha)) + * __temp * std::pow(_Tp{2}, __2alpha) / (__Cp * __C2); + break; + } + if (__its > _S_maxit) + std::__throw_logic_error("__jacobi_zeros: Too many iterations"); + } + __pt[__i - 1].__zero = __z; + __pt[__i - 1].__weight = __w; + } + + return __pt; } } // namespace __detail |