From dd5ed76a707f90f233620805177f6b41b748732c Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@Germany>
Date: Fri, 8 Oct 2010 15:56:58 +0200
Subject: [PATCH] ENH: avoid costly pow() when evaluating Polynomial

    OLD timings:
    ~~~~~~~~~~~~~~~~~~
    evaluate:     -6.82572e+31 in 10.38 s
    hard-coded 0: -6.82572e+31 in 0.2 s
    hard-coded 1: -6.82572e+31 in 10.37 s
    hard-coded 2: -6.82572e+31 in 0.24 s

    New timings:
    ~~~~~~~~~~~~~~~~~~
    evaluate:     -6.82572e+31 in 0.11 s
---
 applications/test/Polynomial/Make/options     |  4 +-
 applications/test/Polynomial/PolynomialTest.C | 94 +++++++++++++++++++
 .../functions/Polynomial/Polynomial.C         |  7 +-
 3 files changed, 101 insertions(+), 4 deletions(-)

diff --git a/applications/test/Polynomial/Make/options b/applications/test/Polynomial/Make/options
index 1b98398577..4c3dd783cb 100644
--- a/applications/test/Polynomial/Make/options
+++ b/applications/test/Polynomial/Make/options
@@ -1,3 +1,3 @@
-EXE_INC = \
+EXE_INC =
 
-EXE_LIBS = \
+EXE_LIBS =
diff --git a/applications/test/Polynomial/PolynomialTest.C b/applications/test/Polynomial/PolynomialTest.C
index 4b7db7ef6f..d3fa03e595 100644
--- a/applications/test/Polynomial/PolynomialTest.C
+++ b/applications/test/Polynomial/PolynomialTest.C
@@ -32,9 +32,14 @@ Description
 #include "IFstream.H"
 #include "Polynomial.H"
 #include "Random.H"
+#include "cpuTime.H"
 
 using namespace Foam;
 
+const int nCoeffs = 8;
+const scalar coeff[] = { 0.11, 0.45, -0.94, 1.58, 2.58, 0.08, 3.15, -4.78 };
+
+
 scalar polyValue(const scalar x)
 {
     // Hard-coded polynomial 8 coeff (7th order)
@@ -65,10 +70,43 @@ scalar intPolyValue(const scalar x)
 }
 
 
+scalar polyValue1(const scalar x)
+{
+    // "normal" evaluation using pow()
+    scalar value = coeff[0];
+
+    for (int i=1; i < nCoeffs; ++i)
+    {
+        value += coeff[i]*pow(x, i);
+    }
+
+    return value;
+}
+
+
+// calculation avoiding pow()
+scalar polyValue2(const scalar x)
+{
+    scalar value = coeff[0];
+
+    scalar powX = x;
+    for (int i=1; i < nCoeffs; ++i, powX *= x)
+    {
+        value += coeff[i] * powX;
+    }
+
+    return value;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
+    const label n = 10000;
+    const label nIters = 1000;
+    scalar sum = 0.0;
+
     IFstream is("polyTestInput");
 
     Polynomial<8> poly("testPoly", is);
@@ -118,6 +156,62 @@ int main(int argc, char *argv[])
         Info<< endl;
     }
 
+
+    //
+    // test speed of Polynomial:
+    //
+    Info<< "start timing loops" << nl
+        << "~~~~~~~~~~~~~~~~~~" << endl;
+
+    cpuTime timer;
+
+    for (int loop = 0; loop < n; ++loop)
+    {
+        sum = 0.0;
+        for (label iter = 0; iter < nIters; ++iter)
+        {
+            sum += poly.evaluate(loop+iter);
+        }
+    }
+    Info<< "evaluate:     " << sum
+        << " in " << timer.cpuTimeIncrement() << " s\n";
+
+
+    for (int loop = 0; loop < n; ++loop)
+    {
+        sum = 0.0;
+        for (label iter = 0; iter < nIters; ++iter)
+        {
+            sum += polyValue(loop+iter);
+        }
+    }
+    Info<< "hard-coded 0: " << sum
+        << " in " << timer.cpuTimeIncrement() << " s\n";
+
+
+    for (int loop = 0; loop < n; ++loop)
+    {
+        sum = 0.0;
+        for (label iter = 0; iter < nIters; ++iter)
+        {
+            sum += polyValue1(loop+iter);
+        }
+    }
+    Info<< "hard-coded 1: " << sum
+        << " in " << timer.cpuTimeIncrement() << " s\n";
+
+    for (int loop = 0; loop < n; ++loop)
+    {
+        sum = 0.0;
+        for (label iter = 0; iter < nIters; ++iter)
+        {
+            sum += polyValue2(loop+iter);
+        }
+    }
+    Info<< "hard-coded 2: " << sum
+        << " in " << timer.cpuTimeIncrement() << " s\n";
+
+
     Info<< nl << "Done." << endl;
 
     return 0;
diff --git a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C
index 1e7b7aa618..903bdffe2b 100644
--- a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C
+++ b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C
@@ -101,9 +101,12 @@ Foam::scalar Foam::Polynomial<PolySize>::evaluate(const scalar x) const
 {
     scalar y = this->v_[0];
 
-    for (label i=1; i<PolySize; i++)
+    // avoid costly pow() in calculation
+    scalar powX = x;
+    for (label i=1; i<PolySize; ++i)
     {
-        y += this->v_[i]*pow(x, i);
+        y += this->v_[i]*powX;
+        powX *= x;
     }
 
     if (logActive_)