]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make a polynomial function inline to allow templated use
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sun, 11 Oct 2020 20:35:22 +0000 (22:35 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 15 Oct 2020 07:50:02 +0000 (09:50 +0200)
include/deal.II/base/polynomial.h
source/base/polynomial.cc

index 1ca3d933796690fe6237c7e812750044aa94f282..1c708c8cd5a2bdc00b6103eeac0c15f89de840eb 100644 (file)
@@ -93,9 +93,11 @@ namespace Polynomials
     /**
      * 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;
@@ -120,14 +122,22 @@ namespace Polynomials
      * 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
@@ -807,6 +817,141 @@ namespace Polynomials
 
 
 
+  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
index 7bf897c29bb65b3374ca0eeeeab2a7c7134543cd..9f9efbeb8de8db7fbe58d3aa59b18339300d77e4 100644 (file)
@@ -106,148 +106,6 @@ namespace Polynomials
 
 
 
-  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()

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.