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

@@ -1,3 +1,3 @@
-EXE_INC = \
+EXE_INC =
 
-EXE_LIBS = \
+EXE_LIBS =

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;

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_)