1 files changed, 6 insertions, 6 deletions
diff --git a/src/nuc_eos_cxx/nuc_eos_press_cs2.cc b/src/nuc_eos_cxx/nuc_eos_press_cs2.cc
index 2652ac5..35de771 100644
--- a/src/nuc_eos_cxx/nuc_eos_press_cs2.cc
+++ b/src/nuc_eos_cxx/nuc_eos_press_cs2.cc
@@ -26,7 +26,7 @@ void CCTK_FNAME(nuc_eos_m_kt1_press_eps_cs2)(const int *restrict n_in,
   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 @@ void CCTK_FNAME(nuc_eos_m_kt1_press_eps_cs2)(const int *restrict n_in,
     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 @@ void CCTK_FNAME(nuc_eos_m_kt0_press_cs2)(const int *restrict n_in,
     // 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 @@ void CCTK_FNAME(nuc_eos_m_kt0_press_cs2)(const int *restrict n_in,
     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 @@ void CCTK_FNAME(nuc_eos_m_kt0_press_cs2)(const int *restrict n_in,
       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 @@ void CCTK_FNAME(nuc_eos_m_kt0_press_cs2)(const int *restrict n_in,
   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
   }