From de5a9ea6c2449409f86907bb287a969eb17064d2 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Fri, 2 Oct 2020 18:03:57 +0200
Subject: [PATCH] Slightly improve evaluation of polynomial

---
 source/base/polynomial.cc | 79 ++++++++++++++++++++++-----------------
 1 file changed, 44 insertions(+), 35 deletions(-)

diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc
index cf1d953307..cbb7bffbdd 100644
--- a/source/base/polynomial.cc
+++ b/source/base/polynomial.cc
@@ -119,6 +119,7 @@ namespace Polynomials
         // 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:
@@ -151,7 +152,7 @@ namespace Polynomials
                 number k_faculty = 1;
                 for (unsigned int k = 0; k <= n_derivatives; ++k)
                   {
-                    values[k] *= k_faculty * lagrange_weight;
+                    values[k] *= k_faculty * weight;
                     k_faculty *= static_cast<number>(k + 1);
                   }
               }
@@ -159,45 +160,53 @@ namespace Polynomials
 
             // 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 loop.
+            // 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:
-              values[0] = 1;
-              for (unsigned int i = 0; i < n_supp; ++i)
-                {
-                  const number v = x - lagrange_support_points[i];
-                  values[0] *= v;
-                }
-              values[0] *= lagrange_weight;
-              break;
+              {
+                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:
-              values[0] = 1;
-              values[1] = 0;
-              for (unsigned int i = 0; i < n_supp; ++i)
-                {
-                  const number v = x - lagrange_support_points[i];
-                  values[1]      = values[1] * v + values[0];
-                  values[0] *= v;
-                }
-              values[0] *= lagrange_weight;
-              values[1] *= lagrange_weight;
-              break;
+              {
+                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:
-              values[0] = 1;
-              values[1] = 0;
-              values[2] = 0;
-              for (unsigned int i = 0; i < n_supp; ++i)
-                {
-                  const number v = x - lagrange_support_points[i];
-                  values[2]      = values[2] * v + values[1];
-                  values[1]      = values[1] * v + values[0];
-                  values[0] *= v;
-                }
-              values[0] *= lagrange_weight;
-              values[1] *= lagrange_weight;
-              values[2] *= static_cast<number>(2) * lagrange_weight;
-              break;
+              {
+                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;
       }
-- 
2.39.5