123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218 |
- #ifndef BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP
- #define BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP
- #include <array>
- #include <cmath>
- #include <limits>
- #include <type_traits>
- namespace boost { namespace math {
- namespace detail {
- template<class Real>
- inline Real B1(Real x)
- {
- if (x < 0)
- {
- return B1(-x);
- }
- if (x < Real(1))
- {
- return 1 - x;
- }
- return Real(0);
- }
- }
- template<unsigned n, typename Real>
- Real cardinal_b_spline(Real x) {
- static_assert(!std::is_integral<Real>::value, "Does not work with integral types.");
- if (x < 0) {
-
- return cardinal_b_spline<n, Real>(-x);
- }
- if (n==0)
- {
- if (x < Real(1)/Real(2)) {
- return Real(1);
- }
- else if (x == Real(1)/Real(2)) {
- return Real(1)/Real(2);
- }
- else {
- return Real(0);
- }
- }
- if (n==1)
- {
- return detail::B1(x);
- }
- Real supp_max = (n+1)/Real(2);
- if (x >= supp_max)
- {
- return Real(0);
- }
-
-
-
- std::array<Real, n> v;
- Real z = x + 1 - supp_max;
- for (unsigned i = 0; i < n; ++i)
- {
- v[i] = detail::B1(z);
- z += 1;
- }
- Real smx = supp_max - x;
- for (unsigned j = 2; j <= n; ++j)
- {
- Real a = (j + 1 - smx);
- Real b = smx;
- for(unsigned k = 0; k <= n - j; ++k)
- {
- v[k] = (a*v[k+1] + b*v[k])/Real(j);
- a += 1;
- b -= 1;
- }
- }
- return v[0];
- }
- template<unsigned n, typename Real>
- Real cardinal_b_spline_prime(Real x)
- {
- static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integer types.");
- if (x < 0)
- {
-
- return -cardinal_b_spline_prime<n, Real>(-x);
- }
- if (n==0)
- {
-
- if (x == Real(1)/Real(2))
- {
- return std::numeric_limits<Real>::infinity();
- }
- else
- {
- return Real(0);
- }
- }
- if (n==1)
- {
- if (x==0)
- {
- return Real(0);
- }
- if (x==1)
- {
- return -Real(1)/Real(2);
- }
- return Real(-1);
- }
- Real supp_max = (n+1)/Real(2);
- if (x >= supp_max)
- {
- return Real(0);
- }
-
- std::array<Real, n> v;
- Real z = x + 1 - supp_max;
- for (unsigned i = 0; i < n; ++i)
- {
- v[i] = detail::B1(z);
- z += 1;
- }
- Real smx = supp_max - x;
- for (unsigned j = 2; j <= n - 1; ++j)
- {
- Real a = (j + 1 - smx);
- Real b = smx;
- for(unsigned k = 0; k <= n - j; ++k)
- {
- v[k] = (a*v[k+1] + b*v[k])/Real(j);
- a += 1;
- b -= 1;
- }
- }
- return v[1] - v[0];
- }
- template<unsigned n, typename Real>
- Real cardinal_b_spline_double_prime(Real x)
- {
- static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integer types.");
- static_assert(n >= 3, "n>=3 for second derivatives of cardinal B-splines is required.");
- if (x < 0)
- {
-
- return cardinal_b_spline_double_prime<n, Real>(-x);
- }
- Real supp_max = (n+1)/Real(2);
- if (x >= supp_max)
- {
- return Real(0);
- }
-
- std::array<Real, n> v;
- Real z = x + 1 - supp_max;
- for (unsigned i = 0; i < n; ++i)
- {
- v[i] = detail::B1(z);
- z += 1;
- }
- Real smx = supp_max - x;
- for (unsigned j = 2; j <= n - 2; ++j)
- {
- Real a = (j + 1 - smx);
- Real b = smx;
- for(unsigned k = 0; k <= n - j; ++k)
- {
- v[k] = (a*v[k+1] + b*v[k])/Real(j);
- a += 1;
- b -= 1;
- }
- }
- return v[2] - 2*v[1] + v[0];
- }
- template<unsigned n, class Real>
- Real forward_cardinal_b_spline(Real x)
- {
- static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integral types.");
- return cardinal_b_spline<n>(x - (n+1)/Real(2));
- }
- }}
- #endif
|