219 lines · plain
1// (C) Copyright Nick Thompson 2019.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP7#define BOOST_MATH_SPECIAL_CARDINAL_B_SPLINE_HPP8 9#include <array>10#include <cmath>11#include <limits>12#include <type_traits>13 14namespace boost { namespace math {15 16namespace 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 33template<unsigned n, typename Real>34Real 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 94template<unsigned n, typename Real>95Real 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 else114 {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 165template<unsigned n, typename Real>166Real 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 210template<unsigned n, class Real>211Real 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#endif219