/**
* Return the value of this polynomial at the given point.
*
- * This function uses the Horner scheme for numerical stability of the
- * evaluation for polynomials in the coefficient form or the product of
- * terms involving the roots if that representation is used.
+ * This function uses the most numerically stable evaluation
+ * algorithm for the provided form of the polynomial. If the
+ * polynomial is in the product form of roots, the evaluation is
+ * based on products of the form (x - x_i), whereas the Horner
+ * scheme is used for polynomials in the coefficient form.
*/
number
value(const number x) const;
* determined by @p n_derivatives and @p values has to provide sufficient
* space for @p n_derivatives + 1 values.
*
- * This function uses the Horner scheme for numerical stability of the
- * evaluation for polynomials in the coefficient form or the product of
- * terms involving the roots if that representation is used.
+ * This function uses the most numerically stable evaluation
+ * algorithm for the provided form of the polynomial. If the
+ * polynomial is in the product form of roots, the evaluation is
+ * based on products of the form (x - x_i), whereas the Horner
+ * scheme is used for polynomials in the coefficient form.
+ *
+ * The template type `Number2` must implement arithmetic
+ * operations such as additions or multiplication with the type
+ * `number` of the polynomial, and must be convertible from
+ * `number` by `operator=`.
*/
+ template <typename Number2>
void
- value(const number x,
+ value(const Number2 x,
const unsigned int n_derivatives,
- number * values) const;
+ Number2 * values) const;
/**
* Degree of the polynomial. This is the degree reflected by the number of
+ template <typename number>
+ template <typename Number2>
+ inline void
+ Polynomial<number>::value(const Number2 x,
+ const unsigned int n_derivatives,
+ Number2 * values) const
+ {
+ // evaluate Lagrange polynomial and derivatives
+ if (in_lagrange_product_form == true)
+ {
+ // to compute the value and all derivatives of a polynomial of the
+ // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
+ // automatic differentiation does.
+ const unsigned int n_supp = lagrange_support_points.size();
+ const number weight = lagrange_weight;
+ switch (n_derivatives)
+ {
+ default:
+ values[0] = 1.;
+ for (unsigned int d = 1; d <= n_derivatives; ++d)
+ values[d] = 0.;
+ for (unsigned int i = 0; i < n_supp; ++i)
+ {
+ const Number2 v = x - lagrange_support_points[i];
+
+ // multiply by (x-x_i) and compute action on all derivatives,
+ // too (inspired from automatic differentiation: implement the
+ // product rule for the old value and the new variable 'v',
+ // i.e., expand value v and derivative one). since we reuse a
+ // value from the next lower derivative from the steps before,
+ // need to start from the highest derivative
+ for (unsigned int k = n_derivatives; k > 0; --k)
+ values[k] = (values[k] * v + values[k - 1]);
+ values[0] *= v;
+ }
+ // finally, multiply by the weight in the Lagrange
+ // denominator. Could be done instead of setting values[0] = 1
+ // above, but that gives different accumulation of round-off
+ // errors (multiplication is not associative) compared to when we
+ // computed the weight, and hence a basis function might not be
+ // exactly one at the center point, which is nice to have. We also
+ // multiply derivatives by k! to transform the product p_n =
+ // p^(n)(x)/k! into the actual form of the derivative
+ {
+ number k_factorial = 1;
+ for (unsigned int k = 0; k <= n_derivatives; ++k)
+ {
+ values[k] *= k_factorial * weight;
+ k_factorial *= static_cast<number>(k + 1);
+ }
+ }
+ break;
+
+ // manually implement case 0 (values only), case 1 (value + first
+ // derivative), and case 2 (up to second derivative) since they
+ // might be called often. then, we can unroll the inner loop and
+ // keep the temporary results as local variables to help the
+ // compiler with the pointer aliasing analysis.
+ case 0:
+ {
+ Number2 value = 1.;
+ for (unsigned int i = 0; i < n_supp; ++i)
+ {
+ const Number2 v = x - lagrange_support_points[i];
+ value *= v;
+ }
+ values[0] = weight * value;
+ break;
+ }
+
+ case 1:
+ {
+ Number2 value = 1.;
+ Number2 derivative = 0.;
+ for (unsigned int i = 0; i < n_supp; ++i)
+ {
+ const Number2 v = x - lagrange_support_points[i];
+ derivative = derivative * v + value;
+ value *= v;
+ }
+ values[0] = weight * value;
+ values[1] = weight * derivative;
+ break;
+ }
+
+ case 2:
+ {
+ Number2 value = 1.;
+ Number2 derivative = 0.;
+ Number2 second = 0.;
+ for (unsigned int i = 0; i < n_supp; ++i)
+ {
+ const Number2 v = x - lagrange_support_points[i];
+ second = second * v + derivative;
+ derivative = derivative * v + value;
+ value *= v;
+ }
+ values[0] = weight * value;
+ values[1] = weight * derivative;
+ values[2] = static_cast<number>(2) * weight * second;
+ break;
+ }
+ }
+ return;
+ }
+
+ Assert(coefficients.size() > 0, ExcEmptyObject());
+
+ // if derivatives are needed, then do it properly by the full
+ // Horner scheme
+ const unsigned int m = coefficients.size();
+ std::vector<Number2> a(coefficients.size());
+ std::copy(coefficients.begin(), coefficients.end(), a.begin());
+ unsigned int j_factorial = 1;
+
+ // loop over all requested derivatives. note that derivatives @p{j>m} are
+ // necessarily zero, as they differentiate the polynomial more often than
+ // the highest power is
+ const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
+ for (unsigned int j = 0; j < min_valuessize_m; ++j)
+ {
+ for (int k = m - 2; k >= static_cast<int>(j); --k)
+ a[k] += x * a[k + 1];
+ values[j] = static_cast<number>(j_factorial) * a[j];
+
+ j_factorial *= j + 1;
+ }
+
+ // fill higher derivatives by zero
+ for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
+ values[j] = 0.;
+ }
+
+
+
template <typename number>
template <class Archive>
inline void
- template <typename number>
- void
- Polynomial<number>::value(const number x,
- const unsigned int n_derivatives,
- number * values) const
- {
- // evaluate Lagrange polynomial and derivatives
- if (in_lagrange_product_form == true)
- {
- // to compute the value and all derivatives of a polynomial of the
- // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
- // automatic differentiation does.
- const unsigned int n_supp = lagrange_support_points.size();
- const number weight = lagrange_weight;
- switch (n_derivatives)
- {
- default:
- values[0] = 1;
- for (unsigned int d = 1; d <= n_derivatives; ++d)
- values[d] = 0;
- for (unsigned int i = 0; i < n_supp; ++i)
- {
- const number v = x - lagrange_support_points[i];
-
- // multiply by (x-x_i) and compute action on all derivatives,
- // too (inspired from automatic differentiation: implement the
- // product rule for the old value and the new variable 'v',
- // i.e., expand value v and derivative one). since we reuse a
- // value from the next lower derivative from the steps before,
- // need to start from the highest derivative
- for (unsigned int k = n_derivatives; k > 0; --k)
- values[k] = (values[k] * v + values[k - 1]);
- values[0] *= v;
- }
- // finally, multiply by the weight in the Lagrange
- // denominator. Could be done instead of setting values[0] = 1
- // above, but that gives different accumulation of round-off
- // errors (multiplication is not associative) compared to when we
- // computed the weight, and hence a basis function might not be
- // exactly one at the center point, which is nice to have. We also
- // multiply derivatives by k! to transform the product p_n =
- // p^(n)(x)/k! into the actual form of the derivative
- {
- number k_factorial = 1;
- for (unsigned int k = 0; k <= n_derivatives; ++k)
- {
- values[k] *= k_factorial * weight;
- k_factorial *= static_cast<number>(k + 1);
- }
- }
- break;
-
- // manually implement case 0 (values only), case 1 (value + first
- // derivative), and case 2 (up to second derivative) since they
- // might be called often. then, we can unroll the inner loop and
- // keep the temporary results as local variables to help the
- // compiler with the pointer aliasing analysis.
- case 0:
- {
- number value = 1;
- for (unsigned int i = 0; i < n_supp; ++i)
- {
- const number v = x - lagrange_support_points[i];
- value *= v;
- }
- values[0] = weight * value;
- break;
- }
-
- case 1:
- {
- number value = 1;
- number derivative = 0;
- for (unsigned int i = 0; i < n_supp; ++i)
- {
- const number v = x - lagrange_support_points[i];
- derivative = derivative * v + value;
- value *= v;
- }
- values[0] = weight * value;
- values[1] = weight * derivative;
- break;
- }
-
- case 2:
- {
- number value = 1;
- number derivative = 0;
- number second = 0;
- for (unsigned int i = 0; i < n_supp; ++i)
- {
- const number v = x - lagrange_support_points[i];
- second = second * v + derivative;
- derivative = derivative * v + value;
- value *= v;
- }
- values[0] = weight * value;
- values[1] = weight * derivative;
- values[2] = static_cast<number>(2) * weight * second;
- break;
- }
- }
- return;
- }
-
- Assert(coefficients.size() > 0, ExcEmptyObject());
-
- // if we only need the value, then call the other function since that is
- // significantly faster (there is no need to allocate and free memory,
- // which is really expensive compared to all the other operations!)
- if (n_derivatives == 0)
- {
- values[0] = value(x);
- return;
- }
-
- // if there are derivatives needed, then do it properly by the full Horner
- // scheme
- const unsigned int m = coefficients.size();
- std::vector<number> a(coefficients);
- unsigned int j_factorial = 1;
-
- // loop over all requested derivatives. note that derivatives @p{j>m} are
- // necessarily zero, as they differentiate the polynomial more often than
- // the highest power is
- const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
- for (unsigned int j = 0; j < min_valuessize_m; ++j)
- {
- for (int k = m - 2; k >= static_cast<int>(j); --k)
- a[k] += x * a[k + 1];
- values[j] = static_cast<number>(j_factorial) * a[j];
-
- j_factorial *= j + 1;
- }
-
- // fill higher derivatives by zero
- for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
- values[j] = 0;
- }
-
-
-
template <typename number>
void
Polynomial<number>::transform_into_standard_form()
projects / dealii.git / commitdiff