boost/math/special_functions/cardinal_b_spline.hpp

   1 //  (C) Copyright Nick Thompson 2019.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP
   7 #define BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP
   8
   9 #include <array>
  10 #include <cmath>
  11 #include <limits>
  12 #include <type_traits>
  13
  14 namespace boost { namespace math {
  15
  16 namespace detail {
  17
  18   template<class Real>
  19   inline Real B1(Real x)
  20   {
  21     if (x < 0)
  22     {
  23       return B1(-x);
  24     }
  25     if (x < Real(1))
  26     {
  27       return 1 - x;
  28     }
  29     return Real(0);
  30   }
  31 }
  32
  33 template<unsigned n, typename Real>
  34 Real cardinal_b_spline(Real x) {
  35     static_assert(!std::is_integral<Real>::value, "Does not work with integral types.");
  36
  37     if (x < 0) {
  38         // All B-splines are even functions:
  39         return cardinal_b_spline<n, Real>(-x);
  40     }
  41
  42     if  (n==0)
  43     {
  44         if (x < Real(1)/Real(2)) {
  45             return Real(1);
  46         }
  47         else if (x == Real(1)/Real(2)) {
  48             return Real(1)/Real(2);
  49         }
  50         else {
  51             return Real(0);
  52         }
  53     }
  54
  55     if (n==1)
  56     {
  57         return detail::B1(x);
  58     }
  59
  60     Real supp_max = (n+1)/Real(2);
  61     if (x >= supp_max)
  62     {
  63         return Real(0);
  64     }
  65
  66     // Fill v with values of B1:
  67     // At most two of these terms are nonzero, and at least 1.
  68     // There is only one non-zero term when n is odd and x = 0.
  69     std::array<Real, n> v;
  70     Real z = x + 1 - supp_max;
  71     for (unsigned i = 0; i < n; ++i)
  72     {
  73         v[i] = detail::B1(z);
  74         z += 1;
  75     }
  76
  77     Real smx = supp_max - x;
  78     for (unsigned j = 2; j <= n; ++j)
  79     {
  80         Real a = (j + 1 - smx);
  81         Real b = smx;
  82         for(unsigned k = 0; k <= n - j; ++k)
  83         {
  84             v[k] = (a*v[k+1] + b*v[k])/Real(j);
  85             a += 1;
  86             b -= 1;
  87         }
  88     }
  89
  90     return v[0];
  91 }
  92
  93
  94 template<unsigned n, typename Real>
  95 Real cardinal_b_spline_prime(Real x)
  96 {
  97     static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integer types.");
  98
  99     if (x < 0)
 100     {
 101         // All B-splines are even functions, so derivatives are odd:
 102         return -cardinal_b_spline_prime<n, Real>(-x);
 103     }
 104
 105
 106     if (n==0)
 107     {
 108         // Kinda crazy but you get what you ask for!
 109         if (x == Real(1)/Real(2))
 110         {
 111             return std::numeric_limits<Real>::infinity();
 112         }
 113         else
 114         {
 115             return Real(0);
 116         }
 117     }
 118
 119     if (n==1)
 120     {
 121         if (x==0)
 122         {
 123             return Real(0);
 124         }
 125         if (x==1)
 126         {
 127             return -Real(1)/Real(2);
 128         }
 129         return Real(-1);
 130     }
 131
 132
 133     Real supp_max = (n+1)/Real(2);
 134     if (x >= supp_max)
 135     {
 136         return Real(0);
 137     }
 138
 139     // 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):
 140     std::array<Real, n> v;
 141     Real z = x + 1 - supp_max;
 142     for (unsigned i = 0; i < n; ++i)
 143     {
 144         v[i] = detail::B1(z);
 145         z += 1;
 146     }
 147
 148     Real smx = supp_max - x;
 149     for (unsigned j = 2; j <= n - 1; ++j)
 150     {
 151         Real a = (j + 1 - smx);
 152         Real b = smx;
 153         for(unsigned k = 0; k <= n - j; ++k)
 154         {
 155             v[k] = (a*v[k+1] + b*v[k])/Real(j);
 156             a += 1;
 157             b -= 1;
 158         }
 159     }
 160
 161     return v[1] - v[0];
 162 }
 163
 164
 165 template<unsigned n, typename Real>
 166 Real cardinal_b_spline_double_prime(Real x)
 167 {
 168     static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integer types.");
 169     static_assert(n >= 3, "n>=3 for second derivatives of cardinal B-splines is required.");
 170
 171     if (x < 0)
 172     {
 173         // All B-splines are even functions, so second derivatives are even:
 174         return cardinal_b_spline_double_prime<n, Real>(-x);
 175     }
 176
 177
 178     Real supp_max = (n+1)/Real(2);
 179     if (x >= supp_max)
 180     {
 181         return Real(0);
 182     }
 183
 184     // 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):
 185     std::array<Real, n> v;
 186     Real z = x + 1 - supp_max;
 187     for (unsigned i = 0; i < n; ++i)
 188     {
 189         v[i] = detail::B1(z);
 190         z += 1;
 191     }
 192
 193     Real smx = supp_max - x;
 194     for (unsigned j = 2; j <= n - 2; ++j)
 195     {
 196         Real a = (j + 1 - smx);
 197         Real b = smx;
 198         for(unsigned k = 0; k <= n - j; ++k)
 199         {
 200             v[k] = (a*v[k+1] + b*v[k])/Real(j);
 201             a += 1;
 202             b -= 1;
 203         }
 204     }
 205
 206     return v[2] - 2*v[1] + v[0];
 207 }
 208
 209
 210 template<unsigned n, class Real>
 211 Real forward_cardinal_b_spline(Real x)
 212 {
 213     static_assert(!std::is_integral<Real>::value, "Cardinal B-splines do not work with integral types.");
 214     return cardinal_b_spline<n>(x - (n+1)/Real(2));
 215 }
 216
 217 }}
 218 #endif