math/special_functions/cardinal_b_spline.hpp

0001 //  (C) Copyright Nick Thompson 2019.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP
0007 #define BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP
0008
0009 #include <array>
0010 #include <cmath>
0011 #include <limits>
0012 #include <type_traits>
0013
0014 namespace boost { namespace math {
0015
0016 namespace detail {
0017
0018   template<class Real>
0019   inline Real B1(Real x)
0020   {
0021     if (x < 0)
0022     {
0023       return B1(-x);
0024     }
0025     if (x < Real(1))
0026     {
0027       return 1 - x;
0028     }
0029     return Real(0);
0030   }
0031 }
0032
0033 template<unsigned n, typename Real>
0034 Real cardinal_b_spline(Real x) {
0035     static_assert(!std::is_integral<Real>::value, "Does not work with integral types.");
0036
0037     if (x < 0) {
0038         // All B-splines are even functions:
0039         return cardinal_b_spline<n, Real>(-x);
0040     }
0041
0042     if  (n==0)
0043     {
0044         if (x < Real(1)/Real(2)) {
0045             return Real(1);
0046         }
0047         else if (x == Real(1)/Real(2)) {
0048             return Real(1)/Real(2);
0049         }
0050         else {
0051             return Real(0);
0052         }
0053     }
0054
0055     if (n==1)
0056     {
0057         return detail::B1(x);
0058     }
0059
0060     Real supp_max = (n+1)/Real(2);
0061     if (x >= supp_max)
0062     {
0063         return Real(0);
0064     }
0065
0066     // Fill v with values of B1:
0067     // At most two of these terms are nonzero, and at least 1.
0068     // There is only one non-zero term when n is odd and x = 0.
0069     std::array<Real, n> v;
0070     Real z = x + 1 - supp_max;
0071     for (unsigned i = 0; i < n; ++i)
0072     {
0073         v[i] = detail::B1(z);
0074         z += 1;
0075     }
0076
0077     Real smx = supp_max - x;
0078     for (unsigned j = 2; j <= n; ++j)
0079     {
0080         Real a = (j + 1 - smx);
0081         Real b = smx;
0082         for(unsigned k = 0; k <= n - j; ++k)
0083         {
0084             v[k] = (a*v[k+1] + b*v[k])/Real(j);
0085             a += 1;
0086             b -= 1;
0087         }
0088     }
0089
0090     return v[0];
0091 }
0092
0093
0094 template<unsigned n, typename Real>
0095 Real cardinal_b_spline_prime(Real x)
0096 {
0097     static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integer types.");
0098
0099     if (x < 0)
0100     {
0101         // All B-splines are even functions, so derivatives are odd:
0102         return -cardinal_b_spline_prime<n, Real>(-x);
0103     }
0104
0105
0106     if (n==0)
0107     {
0108         // Kinda crazy but you get what you ask for!
0109         if (x == Real(1)/Real(2))
0110         {
0111             return std::numeric_limits<Real>::infinity();
0112         }
0113         else
0114         {
0115             return Real(0);
0116         }
0117     }
0118
0119     if (n==1)
0120     {
0121         if (x==0)
0122         {
0123             return Real(0);
0124         }
0125         if (x==1)
0126         {
0127             return -Real(1)/Real(2);
0128         }
0129         return Real(-1);
0130     }
0131
0132
0133     Real supp_max = (n+1)/Real(2);
0134     if (x >= supp_max)
0135     {
0136         return Real(0);
0137     }
0138
0139     // Now we want to evaluate B_{n}(x), but stop at the second to last step and collect B_{n-1}(x+1/2) and B_{n-1}(x-1/2):
0140     std::array<Real, n> v;
0141     Real z = x + 1 - supp_max;
0142     for (unsigned i = 0; i < n; ++i)
0143     {
0144         v[i] = detail::B1(z);
0145         z += 1;
0146     }
0147
0148     Real smx = supp_max - x;
0149     for (unsigned j = 2; j <= n - 1; ++j)
0150     {
0151         Real a = (j + 1 - smx);
0152         Real b = smx;
0153         for(unsigned k = 0; k <= n - j; ++k)
0154         {
0155             v[k] = (a*v[k+1] + b*v[k])/Real(j);
0156             a += 1;
0157             b -= 1;
0158         }
0159     }
0160
0161     return v[1] - v[0];
0162 }
0163
0164
0165 template<unsigned n, typename Real>
0166 Real cardinal_b_spline_double_prime(Real x)
0167 {
0168     static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integer types.");
0169     static_assert(n >= 3, "n>=3 for second derivatives of cardinal B-splines is required.");
0170
0171     if (x < 0)
0172     {
0173         // All B-splines are even functions, so second derivatives are even:
0174         return cardinal_b_spline_double_prime<n, Real>(-x);
0175     }
0176
0177
0178     Real supp_max = (n+1)/Real(2);
0179     if (x >= supp_max)
0180     {
0181         return Real(0);
0182     }
0183
0184     // Now we want to evaluate B_{n}(x), but stop at the second to last step and collect B_{n-1}(x+1/2) and B_{n-1}(x-1/2):
0185     std::array<Real, n> v;
0186     Real z = x + 1 - supp_max;
0187     for (unsigned i = 0; i < n; ++i)
0188     {
0189         v[i] = detail::B1(z);
0190         z += 1;
0191     }
0192
0193     Real smx = supp_max - x;
0194     for (unsigned j = 2; j <= n - 2; ++j)
0195     {
0196         Real a = (j + 1 - smx);
0197         Real b = smx;
0198         for(unsigned k = 0; k <= n - j; ++k)
0199         {
0200             v[k] = (a*v[k+1] + b*v[k])/Real(j);
0201             a += 1;
0202             b -= 1;
0203         }
0204     }
0205
0206     return v[2] - 2*v[1] + v[0];
0207 }
0208
0209
0210 template<unsigned n, class Real>
0211 Real forward_cardinal_b_spline(Real x)
0212 {
0213     static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integral types.");
0214     return cardinal_b_spline<n>(x - (n+1)/Real(2));
0215 }
0216
0217 }}
0218 #endif