Loading [MathJax]/extensions/tex2jax.js
LSST Applications g0fba68d861+83433b07ee,g16d25e1f1b+23bc9e47ac,g1ec0fe41b4+3ea9d11450,g1fd858c14a+9be2b0f3b9,g2440f9efcc+8c5ae1fdc5,g35bb328faa+8c5ae1fdc5,g4a4af6cd76+d25431c27e,g4d2262a081+c74e83464e,g53246c7159+8c5ae1fdc5,g55585698de+1e04e59700,g56a49b3a55+92a7603e7a,g60b5630c4e+1e04e59700,g67b6fd64d1+3fc8cb0b9e,g78460c75b0+7e33a9eb6d,g786e29fd12+668abc6043,g8352419a5c+8c5ae1fdc5,g8852436030+60e38ee5ff,g89139ef638+3fc8cb0b9e,g94187f82dc+1e04e59700,g989de1cb63+3fc8cb0b9e,g9d31334357+1e04e59700,g9f33ca652e+0a83e03614,gabe3b4be73+8856018cbb,gabf8522325+977d9fabaf,gb1101e3267+8b4b9c8ed7,gb89ab40317+3fc8cb0b9e,gc0af124501+57ccba3ad1,gcf25f946ba+60e38ee5ff,gd6cbbdb0b4+1cc2750d2e,gd794735e4e+7be992507c,gdb1c4ca869+be65c9c1d7,gde0f65d7ad+c7f52e58fe,ge278dab8ac+6b863515ed,ge410e46f29+3fc8cb0b9e,gf35d7ec915+97dd712d81,gf5e32f922b+8c5ae1fdc5,gf618743f1b+747388abfa,gf67bdafdda+3fc8cb0b9e,w.2025.18
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
RecurrenceBasis1d.h
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2/*
3 * Developed for the LSST Data Management System.
4 * This product includes software developed by the LSST Project
5 * (https://www.lsst.org).
6 * See the COPYRIGHT file at the top-level directory of this distribution
7 * for details of code ownership.
8 *
9 * This program is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 3 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program. If not, see <https://www.gnu.org/licenses/>.
21 */
22#ifndef LSST_AFW_MATH_POLYNOMIALS_RecurrenceBasis1d_h_INCLUDED
23#define LSST_AFW_MATH_POLYNOMIALS_RecurrenceBasis1d_h_INCLUDED
24
27
28namespace lsst { namespace geom { namespace polynomials {
29
30template <typename Basis>
31class Function1d;
32
33#ifdef DOXYGEN
34
43struct Recurrence {
44
46 static double getB0(double x);
47
49 static double getB1(double x);
50
61 static double next(double x, std::size_t n, double current, double previous);
62
63};
64
65#endif // DOXYGEN
66
67
84template <typename Recurrence>
86public:
87
90
93
95 explicit RecurrenceBasis1d(std::size_t order) noexcept :
96 _order(order)
97 {}
98
101
104
107
110
112 std::size_t getOrder() const noexcept { return _order; }
113
115 std::size_t size() const noexcept { return _order + 1; }
116
123 Scaled scaled(Scaling1d const & scaling) const noexcept {
124 return Scaled(*this, scaling);
125 }
126
145 template <typename Vector>
146 double sumWith(double x, Vector const & coefficients, SumMode mode=SumMode::FAST) const {
147 // This universal lambda lets us effectively template most of the
148 // implementation of this function on double vs. SafeSum<double>
149 // without having to define an external template.
150 auto accumulate = [x, coefficients, this](auto & sum) {
151 double previous = Recurrence::getB0(x);
152 if (_order > 0u) {
153 double current = Recurrence::getB1(x);
154 sum += coefficients[1]*current;
155 for (std::size_t n = 2; n <= _order; ++n) {
156 double next = Recurrence::next(x, n, current, previous);
157 sum += coefficients[n]*next;
158 previous = current;
159 current = next;
160 }
161 }
162 };
163 double result = 0.0;
164 if (mode == SumMode::FAST) {
165 double z = Recurrence::getB0(x)*coefficients[0];
166 accumulate(z);
167 result = z;
168 } else {
169 SafeSum<double> z(Recurrence::getB0(x)*coefficients[0]);
170 accumulate(z);
171 result = static_cast<double>(z);
172 }
173 return result;
174 }
175
186 template <typename Vector>
187 void fill(double x, Vector && basis) const {
189 if (_order > 0u) {
191 for (std::size_t n = 2; n <= _order; ++n) {
193 x, n,
194 std::forward<Vector>(basis)[n - 1],
195 std::forward<Vector>(basis)[n - 2]
196 );
197 }
198 }
199 }
200
201private:
202 std::size_t _order;
203};
204
205}}} // namespace lsst::geom::polynomials
206
207#endif // !LSST_AFW_MATH_POLYNOMIALS_RecurrenceBasis1d_h_INCLUDED
A 1-d function defined by a series expansion and its coefficients.
Definition Function1d.h:42
RecurrenceBasis1d & operator=(RecurrenceBasis1d &&)=default
Default move assignment.
RecurrenceBasis1d(RecurrenceBasis1d const &)=default
Default copy constructor.
RecurrenceBasis1d & operator=(RecurrenceBasis1d const &)=default
Default copy assignment.
RecurrenceBasis1d(std::size_t order) noexcept
Construct a basis with the given order (inclusive).
void fill(double x, Vector &&basis) const
Evaluate the basis at a given point.
Scaled scaled(Scaling1d const &scaling) const noexcept
Return a scaled basis with the same order and recurrence.
double sumWith(double x, Vector const &coefficients, SumMode mode=SumMode::FAST) const
Evaluate a basis expansion with the given coefficients.
RecurrenceBasis1d(RecurrenceBasis1d &&)=default
Default move constructor.
std::size_t getOrder() const noexcept
Return the order of the basis.
std::size_t size() const noexcept
Return the number of elements in the basis.
A numerically stable summation algorithm for floating-point numbers.
Definition SafeSum.h:62
A 1-d basis that transforms all input points before evaluating nested basis.
A 1-d affine transform that can be used to map one interval to another.
Definition Scaling1d.h:46
T forward(T... args)
Low-level polynomials (including special polynomials) in C++.
SumMode
Enum used to control how to sum polynomial terms.
Definition SafeSum.h:32
@ FAST
Summation using regular floating-point addition.
Definition SafeSum.h:34
A recurrence relation concept for RecurrenceBasis1d.
static double getB1(double x)
Return the first element of the basis, .
static double next(double x, std::size_t n, double current, double previous)
Return the next element in the recurrence.
static double getB0(double x)
Return the zeroth element of the basis, .