[svn:einsteintoolkit] EOS_Omni/trunk/src/nuc_eos_cxx/ (Rev. 101) - Commits

12 Mar 2014

User: rhaas
Date: 2014/03/12 10:01 PM
Modified:
 /trunk/src/nuc_eos_cxx/
  helpers.hh, nuc_eos_dpdrhoe_dpderho.cc, nuc_eos_full.cc, nuc_eos_press.cc, nuc_eos_press_cs2.cc, nuc_eos_short.cc, readtable.cc
Log:
 EOS_Omni: low-level optimizations
* use log-rules to transform base 10 logs to natural ones
   faster and hopefully more accurate
 * optimize away a division when computing cs2
 * optimize bracketing test to (a-a1)*(a-a2)<0
 * replace divisions by multiplication by inverse
 * use CCTK_BUILTIN_EXPECT to indicate likely outcome
   this should help branch prediction since the compiler can put the the
   less likely branch into the unfavorable location.
 * remove superfluous if statements
   whenever we leave a itmax loop regularly we iterated for too long
 * explicitly replace x/y by x*1./y
From: Roland Haas rhaas@tapir.caltech.edu
File Changes:
Directory: /trunk/src/nuc_eos_cxx/
==================================
File [modified]: helpers.hh
Delta lines: +49 -52
===================================================================

--- trunk/src/nuc_eos_cxx/helpers.hh	2014-03-13 03:01:52 UTC (rev 100)
+++ trunk/src/nuc_eos_cxx/helpers.hh	2014-03-13 03:01:55 UTC (rev 101)
@@ -18,24 +18,24 @@
   // 105 -- rho too high
   // 106 -- rho too low
-  if(xrho > eos_rhomax) {
+  if(CCTK_BUILTIN_EXPECT(xrho > eos_rhomax, false)) {
     return 105;
   }
-  if(xrho < eos_rhomin) {
+  if(CCTK_BUILTIN_EXPECT(xrho < eos_rhomin, false)) {
     return 106;
   }
-  if(xye > eos_yemax) {
+  if(CCTK_BUILTIN_EXPECT(xye > eos_yemax, false)) {
     return 101;
   }
-  if(xye < eos_yemin) {
+  if(CCTK_BUILTIN_EXPECT(xye < eos_yemin, false)) {
     // this is probably not pure and should be removed
     fprintf(stderr,"xye: %15.6E eos_yemin: %15.6E\n",xye,eos_yemin);
     return 102;
   }
-  if(xtemp > eos_tempmax) {
+  if(CCTK_BUILTIN_EXPECT(xtemp > eos_tempmax, false)) {
     return 103;
   }
-  if(xtemp < eos_tempmin) {
+  if(CCTK_BUILTIN_EXPECT(xtemp < eos_tempmin, false)) {
     return 104;
   }
   return 0;
@@ -53,16 +53,16 @@
   // 105 -- rho too high
   // 106 -- rho too low
-  if(xrho > eos_rhomax) {
+  if(CCTK_BUILTIN_EXPECT(xrho > eos_rhomax, false)) {
     return 105;
   }
-  if(xrho < eos_rhomin) {
+  if(CCTK_BUILTIN_EXPECT(xrho < eos_rhomin, false)) {
     return 106;
   }
-  if(xye > eos_yemax) {
+  if(CCTK_BUILTIN_EXPECT(xye > eos_yemax, false)) {
     return 101;
   }
-  if(xye < eos_yemin) {
+  if(CCTK_BUILTIN_EXPECT(xye < eos_yemin, false)) {
     return 102;
   }
   return 0;
@@ -331,10 +331,12 @@
   // then interpolate in y
   // assume rectangular grid
-  double t1 = (fs[1]-fs[0])/(xs[1]-xs[0]) * (x - xs[0]) + fs[0];
-  double t2 = (fs[3]-fs[2])/(xs[1]-xs[0]) * (x - xs[0]) + fs[2];
+  double dxi = 1./(xs[1]-xs[0]);
+  double dyi = 1./(ys[1]-ys[0]); // x*1./y uses faster instructions than x/y
+  double t1 = (fs[1]-fs[0])*dxi * (x - xs[0]) + fs[0];
+  double t2 = (fs[3]-fs[2])*dxi * (x - xs[0]) + fs[2];
-  return (t2 - t1)/(ys[1]-ys[0]) * (y-ys[0]) + t1;
+  return (t2 - t1)*dyi * (y-ys[0]) + t1;
 }
static inline __attribute__((always_inline))
@@ -359,6 +361,8 @@
   const double dltp = log(1.2);
   const double dltm = log(0.8);
+  const double leps0_prec = leps0*prec;
+
   // temporary local vars
   double lt, lt1, lt2;
   double ltmin = logtemp[0];
@@ -404,7 +408,7 @@
 #endif
bcount++;
-    if(bcount >= maxbcount) {
+    if(CCTK_BUILTIN_EXPECT(bcount >= maxbcount, false)) {
       *keyerrt = 667;
       return;
     }
@@ -432,22 +436,20 @@
     nuc_eos_C_linterp_one(idx,delx,dely,delz,&f2a,iv);
fmid=f2a-leps0;
-    if(fmid <= 0.0) lt=ltmid;
+    if(CCTK_BUILTIN_EXPECT(fmid <= 0.0, false)) lt=ltmid;
 #if DEBUG
     fprintf(stderr,"bisection step 2 it %d, fmid: %15.6E f2a: %15.6E lt: %15.6E ltmid: %15.6E dlt: %15.6E\n",
        it,fmid,f2a,lt,ltmid,dlt);
 #endif
-    if(fabs(1.0-f2a/leps0) <= prec) {
+    if(CCTK_BUILTIN_EXPECT(fabs(leps0-f2a) <= leps0_prec, false)) {
       *ltout = ltmid;
       return;
     }
   } // for it = 0
-  if(it >= itmax-1) {
-    *keyerrt = 667;
-    return;
-  }
+  *keyerrt = 667;
+  return;
 } // bisection
@@ -465,7 +467,7 @@
// local variables
   const int itmax = 200; // use at most 10 iterations, then go to bisection
-  double dlepsdlt; // derivative dlogeps/dlogT
+  double dlepsdlti; // 1 / derivative dlogeps/dlogT
   double ldt;
   double leps,leps0; // temp vars for eps
   double ltn, lt; // temp vars for temperature
@@ -489,6 +491,7 @@
   fprintf(stderr,"it: %d t: %15.6E leps: %15.6E eps0: %15.6E del: %15.6E\n",
      it,lt, leps,leps0,fabs(leps-leps0)/(fabs(leps0)));
 #endif
+  // TODO: profile this to see which outcome is more likely
   if(fabs(leps-leps0) < prec*fabs(leps0)) {
     *ltout = lt0;
     return;
@@ -548,8 +551,7 @@
// Check if we are already bracketing the input internal
     // energy. If so, interpolate for new T.
-    if( (leps0 >= epst1 && leps0 <= epst2) 
-	|| (leps0 <= epst1 && leps0 >=epst2)) {
+    if(CCTK_BUILTIN_EXPECT((leps0 - epst1) * (leps0 - epst2) <= 0., false)) {
*ltout = (logtemp[itemp]-logtemp[itemp-1]) / (epst2 - epst1) * 
    (leps0 - epst1) + logtemp[itemp-1];
@@ -559,8 +561,8 @@
// well, then do a Newton-Raphson step
     // first, guess the derivative
-    dlepsdlt = (epst2-epst1)/(logtemp[itemp]-logtemp[itemp-1]);
-    ldt = -(leps - leps0) / dlepsdlt * fac;
+    dlepsdlti = (logtemp[itemp]-logtemp[itemp-1])/(epst2-epst1);
+    ldt = -(leps - leps0) * dlepsdlti * fac;
ltn = MIN(MAX(lt + ldt,ltmin),ltmax);
     lt = ltn;
@@ -578,7 +580,7 @@
     if(oerr < err) fac *= 0.9;
     oerr = err;
-    if(err < prec*fabs(leps0)) {
+    if(CCTK_BUILTIN_EXPECT(err < prec*fabs(leps0), false)) {
       *ltout = lt;
       return;
     }
@@ -587,20 +589,18 @@
} // while(it < itmax)
-  if(it >= itmax - 1) {
-    // try bisection
-    // bisection(lr, lt0, ye, leps0, ltout, 1, prec, keyerrt);
+  // try bisection
+  // bisection(lr, lt0, ye, leps0, ltout, 1, prec, keyerrt);
 #if DEBUG
-    fprintf(stderr, "Failed to converge. This is bad. Trying bisection!\n");
+  fprintf(stderr, "Failed to converge. This is bad. Trying bisection!\n");
 #endif
-    bisection(lr,lt0,ye,leps0,prec,ltout,1,keyerrt);
+  bisection(lr,lt0,ye,leps0,prec,ltout,1,keyerrt);
 #if DEBUG
-    if(*keyerrt==667) {
-      fprintf(stderr,"This is worse. Bisection failed!\n");
-      abort();
-    }      
+  if(*keyerrt==667) {
+    fprintf(stderr,"This is worse. Bisection failed!\n");
+    abort();
+  }      
 #endif
-  }
return;
@@ -620,7 +620,7 @@
// local variables
     const int itmax = 200; // use at most 10 iterations, then go to bisection
-    double dentdlt; // derivative dentropy/dlogT
+    double dentdlti; // 1 / derivative dentropy/dlogT
     double ldt;
     double ent,ent0; // temp vars for eps
     double ltn, lt; // temp vars for temperature
@@ -644,7 +644,7 @@
     fprintf(stderr,"it: %d t: %15.6E leps: %15.6E eps0: %15.6E del: %15.6E\n",
        it,lt,ent,ent0,fabs(ent-ent0)/(fabs(ent0)));
 #endif
-    if(fabs(ent-ent0) < prec*fabs(ent0)) {
+    if(CCTK_BUILTIN_EXPECT(fabs(ent-ent0) < prec*fabs(ent0), false)) {
       *ltout = lt0;
       return;
     }
@@ -703,8 +703,7 @@
// Check if we are already bracketing the input internal
       // energy. If so, interpolate for new T.
-      if( (ent0 >= ent1 && ent0 <= ent2) 
-	  || (ent0 <= ent1 && ent0 >= ent2)) {
+      if(CCTK_BUILTIN_EXPECT((ent0 - ent1) * (ent0 - ent2) <= 0., false)) {
*ltout = (logtemp[itemp]-logtemp[itemp-1]) / (ent2 - ent1) * 
      (ent0 - ent1) + logtemp[itemp-1];
@@ -714,8 +713,8 @@
// well, then do a Newton-Raphson step
       // first, guess the derivative
-      dentdlt = (ent2-ent1)/(logtemp[itemp]-logtemp[itemp-1]);
-      ldt = -(ent - ent0) / dentdlt * fac;
+      dentdlti = (logtemp[itemp]-logtemp[itemp-1])/(ent2-ent1);
+      ldt = -(ent - ent0) * dentdlti * fac;
ltn = MIN(MAX(lt + ldt,ltmin),ltmax);
       lt = ltn;
@@ -733,7 +732,7 @@
       if(oerr < err) fac *= 0.9;
       oerr = err;
-      if(err < prec*fabs(ent0)) {
+      if(CCTK_BUILTIN_EXPECT(err < prec*fabs(ent0), false)) {
    *ltout = lt;
    return;
       }
@@ -742,19 +741,17 @@
} // while(it < itmax)
-    if(it >= itmax - 1) {
-      // try bisection
+    // try bisection
 #if DEBUG
-      fprintf(stderr, "Failed to converge. This is bad. Trying bisection!\n");
+    fprintf(stderr, "Failed to converge. This is bad. Trying bisection!\n");
 #endif
-      bisection(lr,lt0,ye,ent0,prec,ltout,2,keyerrt);
+    bisection(lr,lt0,ye,ent0,prec,ltout,2,keyerrt);
 #if DEBUG
-      if(*keyerrt==667) {
-	fprintf(stderr,"This is worse. Bisection failed!\n");
-	abort();
-      }      
+    if(*keyerrt==667) {
+      fprintf(stderr,"This is worse. Bisection failed!\n");
+      abort();
+    }      
 #endif
-    }
return;
File [modified]: nuc_eos_dpdrhoe_dpderho.cc
Delta lines: +3 -3
===================================================================
--- trunk/src/nuc_eos_cxx/nuc_eos_dpdrhoe_dpderho.cc	2014-03-13 03:01:52 UTC (rev 100)
+++ trunk/src/nuc_eos_cxx/nuc_eos_dpdrhoe_dpderho.cc	2014-03-13 03:01:55 UTC (rev 101)
@@ -30,7 +30,7 @@
     // Note that this code now requires that the
     // temperature guess be within the table bounds
     keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]);
-    if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     }
   }
@@ -45,7 +45,7 @@
     const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax));
     double ltout;
     const double epstot = eps[i]+energy_shift;
-    if(epstot>0.0e0) {
+    if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) {
       // this is the standard scenario; eps is larger than zero
       // and we can operate with logarithmic tables
       const double lxeps = log(epstot);
@@ -61,7 +61,7 @@
       keyerr[i] = 667;
     } // epstot > 0.0
-    if(keyerr[i]) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       // now try negative temperature treatment
       double eps0, eps1;
       int idx[8];
File [modified]: nuc_eos_full.cc
Delta lines: +6 -6
===================================================================
--- trunk/src/nuc_eos_cxx/nuc_eos_full.cc	2014-03-13 03:01:52 UTC (rev 100)
+++ trunk/src/nuc_eos_cxx/nuc_eos_full.cc	2014-03-13 03:01:55 UTC (rev 101)
@@ -39,7 +39,7 @@
   for(int i=0;i<n;i++) {
     // check if we are fine
     keyerr[i] = checkbounds(rho[i], temp[i], ye[i]);
-    if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     }
   }
@@ -148,7 +148,7 @@
     prs[i] = exp(prs[i]);
     eps[i] = exp(eps[i]) - energy_shift;
 #if HAVEGR
-    cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]);
+    cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]);
 #endif
   }
@@ -194,7 +194,7 @@
     // Note that this code now requires that the
     // temperature guess be within the table bounds
     keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]);
-    if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     }
   }
@@ -209,7 +209,7 @@
     const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax));
     double ltout;
     const double epstot = eps[i]+energy_shift;
-    if(epstot>0.0e0) {
+    if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) {
       // this is the standard scenario; eps is larger than zero
       // and we can operate with logarithmic tables
       const double lxeps = log(epstot);
@@ -225,7 +225,7 @@
     } // epstot > 0.0
-    if(keyerr[i]) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr=1;
     } else {
       temp[i] = exp(ltout);
@@ -318,7 +318,7 @@
   for(int i=0;i<n;i++) {
     prs[i] = exp(prs[i]);
 #if HAVEGR
-    cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]);
+    cs2[i] = rho[i]*cs2[i] / (rho[i] + rho[i]*eps[i] + prs[i]);
 #endif
   }
File [modified]: nuc_eos_press.cc
Delta lines: +4 -4
===================================================================
--- trunk/src/nuc_eos_cxx/nuc_eos_press.cc	2014-03-13 03:01:52 UTC (rev 100)
+++ trunk/src/nuc_eos_cxx/nuc_eos_press.cc	2014-03-13 03:01:55 UTC (rev 101)
@@ -25,7 +25,7 @@
   for(int i=0;i<n;i++) {
     // check if we are fine
     keyerr[i] = checkbounds(rho[i], temp[i], ye[i]);
-    if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     }
   }
@@ -87,7 +87,7 @@
     // Note that this code now requires that the
     // temperature guess be within the table bounds
     keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]);
-    if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     }
   }
@@ -103,7 +103,7 @@
     const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax));
     double ltout;
     const double epstot = eps[i]+energy_shift;
-    if(epstot>0.0e0) {
+    if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) {
       // this is the standard scenario; eps is larger than zero
       // and we can operate with logarithmic tables
       const double lxeps = log(epstot);
@@ -119,7 +119,7 @@
       keyerr[i] = 667;
     } // epstot > 0.0
-    if(keyerr[i]) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       // now try negative temperature treatment
       double eps0, eps1;
       int idx[8];
File [modified]: nuc_eos_press_cs2.cc
Delta lines: +6 -6
===================================================================
--- trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc	2014-03-13 03:01:52 UTC (rev 100)
+++ trunk/src/nuc_eos_cxx/nuc_eos_press_cs2.cc	2014-03-13 03:01:55 UTC (rev 101)
@@ -26,7 +26,7 @@
   for(int i=0;i<n;i++) {
     // check if we are fine
     keyerr[i] = checkbounds(rho[i], temp[i], ye[i]);
-    if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     }
   }
@@ -63,7 +63,7 @@
     prs[i] = exp(prs[i]);
     eps[i] = exp(eps[i]) - energy_shift;
 #if HAVEGR
-    cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]);
+    cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]);
 #endif
   }
@@ -96,7 +96,7 @@
     // Note that this code now requires that the
     // temperature guess be within the table bounds
     keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]);
-    if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     }
   }
@@ -112,7 +112,7 @@
     const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax));
     double ltout;
     const double epstot = eps[i]+energy_shift;
-    if(epstot>0.0e0) {
+    if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) {
       // this is the standard scenario; eps is larger than zero
       // and we can operate with logarithmic tables
       const double lxeps = log(epstot);
@@ -124,7 +124,7 @@
       keyerr[i] = 667;
     } // epstot > 0.0
-    if(keyerr[i]) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       // now try negative temperature treatment
       double eps0, eps1;
       int idx[8];
@@ -175,7 +175,7 @@
   for(int i=0;i<n;i++) {
     prs[i] = exp(prs[i]);
 #if HAVEGR
-    cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]);
+    cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i]*eps[i] + prs[i]);
 #endif
   }
File [modified]: nuc_eos_short.cc
Delta lines: +9 -9
===================================================================
--- trunk/src/nuc_eos_cxx/nuc_eos_short.cc	2014-03-13 03:01:52 UTC (rev 100)
+++ trunk/src/nuc_eos_cxx/nuc_eos_short.cc	2014-03-13 03:01:55 UTC (rev 101)
@@ -31,7 +31,7 @@
   for(int i=0;i<n;i++) {
     // check if we are fine
     keyerr[i] = checkbounds(rho[i], temp[i], ye[i]);
-    if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     }
   }
@@ -94,7 +94,7 @@
     prs[i] = exp(prs[i]);
     eps[i] = exp(eps[i]) - energy_shift;
 #if HAVEGR
-    cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]);
+    cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i]*eps[i] + prs[i]);
 #endif
   }
@@ -131,7 +131,7 @@
     // Note that this code now requires that the
     // temperature guess be within the table bounds
     keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]);
-    if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     }
   }
@@ -146,7 +146,7 @@
     const double lt = log(MIN(MAX(temp[i],eos_tempmin),eos_tempmax));
     double ltout;
     const double epstot = eps[i]+energy_shift;
-    if(epstot>0.0e0) {
+    if(CCTK_BUILTIN_EXPECT(epstot>0.0e0, true)) {
       // this is the standard scenario; eps is larger than zero
       // and we can operate with logarithmic tables
       const double lxeps = log(epstot);
@@ -163,7 +163,7 @@
     } // epstot > 0.0
-    if(keyerr[i]) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
       *anyerr = 1;
     } else {
       temp[i] = exp(ltout);
@@ -211,7 +211,7 @@
   for(int i=0;i<n;i++) {
     prs[i] = exp(prs[i]);
 #if HAVEGR
-    cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]);
+    cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]);
 #endif
   }
@@ -247,7 +247,7 @@
       // Note that this code now requires that the
       // temperature guess be within the table bounds
       keyerr[i] = checkbounds_kt0_noTcheck(rho[i], ye[i]);
-      if(keyerr[i] != 0) {
+    if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
    *anyerr = 1;
       }
     }
@@ -264,7 +264,7 @@
       nuc_eos_findtemp_entropy(lr,lt,ye[i],ent[i],*prec,
    		       (double *restrict)(&ltout),&keyerr[i]);
-      if(keyerr[i]) {
+      if(CCTK_BUILTIN_EXPECT(keyerr[i] != 0, false)) {
         *anyerr = 1;
       } else {
    temp[i] = exp(ltout);
@@ -313,7 +313,7 @@
       prs[i] = exp(prs[i]);
       eps[i] = exp(eps[i]) - energy_shift;
 #if HAVEGR
-      cs2[i] = cs2[i] / (1.0 + eps[i] + prs[i]/rho[i]);
+      cs2[i] = rho[i] * cs2[i] / (rho[i] + rho[i] * eps[i] + prs[i]);
 #endif
     }
File [modified]: readtable.cc
Delta lines: +8 -4
===================================================================
--- trunk/src/nuc_eos_cxx/readtable.cc	2014-03-13 03:01:52 UTC (rev 100)
+++ trunk/src/nuc_eos_cxx/readtable.cc	2014-03-13 03:01:55 UTC (rev 101)
@@ -194,11 +194,15 @@
   // pressure
   energy_shift = energy_shift * EPSGF;
   for(int i=0;i<nrho;i++) {
-    logrho[i] = log(pow(10.0,logrho[i]) * RHOGF);
+    // rewrite:
+    //logrho[i] = log(pow(10.0,logrho[i]) * RHOGF);
+    // by using log(a^b*c) = b*log(a)+log(c)
+    logrho[i] = logrho[i] * log(10.) + log(RHOGF);
   }
for(int i=0;i<ntemp;i++) {
-    logtemp[i] = log(pow(10.0,logtemp[i]));
+    //logtemp[i] = log(pow(10.0,logtemp[i]));
+    logtemp[i] = logtemp[i]*log(10.0);
   }
// allocate epstable; a linear-scale eps table
@@ -214,12 +218,12 @@
{ // pressure
       int idx = 0 + NTABLES*i;
-      alltables[idx] = log(pow(10.0,alltables[idx])*PRESSGF);
+      alltables[idx] = alltables[idx] * log(10.0) + log(PRESSGF);
     }
{ // eps
       int idx = 1 + NTABLES*i;
-      alltables[idx] = log(pow(10.0,alltables[idx])*EPSGF);
+      alltables[idx] = alltables[idx] * log(10.0) + log(EPSGF);
       epstable[i] = exp(alltables[idx]);
     }

    