   t = mul_fp64(r, INVERSE_FACTORIAL_5_FP64);
   s = sum_fp64(s, t);
 
-  r = mul_fp64(r, x);
-  t = mul_fp64(r, INVERSE_FACTORIAL_7_FP64);
-  s = sum_fp64(s, t);
-
-  r = mul_fp64(r, x);
-  t = mul_fp64(r, INVERSE_FACTORIAL_9_FP64);
-  s = sum_fp64(s, t);
+  /* keep the following commented code in case we need them
+  for extra accuracy from the Taylor expansion*/
+  // r = mul_fp64(r, x);
+  // t = mul_fp64(r, INVERSE_FACTORIAL_7_FP64);
+  // s = sum_fp64(s, t);
+
+  // r = mul_fp64(r, x);
+  // t = mul_fp64(r, INVERSE_FACTORIAL_9_FP64);
+  // s = sum_fp64(s, t);
 
   return s;
 }

  p = mul_fp64(p, r);
  t = mul_fp64(p, INVERSE_FACTORIAL_5_FP64);

  s = sum_fp64(s, t);
  p = mul_fp64(p, r);
  t = mul_fp64(p, INVERSE_FACTORIAL_6_FP64);

  s = sum_fp64(s, t);
  p = mul_fp64(p, r);
  t = mul_fp64(p, INVERSE_FACTORIAL_7_FP64);

  s = sum_fp64(s, t);

  // At this point, s = exp(r) - 1; but after following 4 recursions, we will get exp(r) ^ 512 - 1.

  s = sum_fp64(s * 2.0, mul_fp64(s, s));
  s = sum_fp64(s * 2.0, mul_fp64(s, s));
  s = sum_fp64(s * 2.0, mul_fp64(s, s));
  s = sum_fp64(s * 2.0, mul_fp64(s, s));

  // We can add more iterations here and increase k.
  s = sum_fp64(s * 2.0, mul_fp64(s, s));
  s = sum_fp64(s * 2.0, mul_fp64(s, s));
  s = sum_fp64(s * 2.0, mul_fp64(s, s));
  s = sum_fp64(s * 2.0, mul_fp64(s, s));
  s = sum_fp64(s * 2.0, mul_fp64(s, s));

#if defined(NVIDIA_FP64_WORKAROUND) || defined(INTEL_FP64_WORKAROUND)
  s = sum_fp64(s, vec2(1.0, 0.0) * ONE);
#else
  s = sum_fp64(s, vec2(1.0, 0.0));
#endif

#endif

#if defined(NVIDIA_EQUATION_WORKAROUND) || defined(INTEL_EQUATION_WORKAROUND)
/* Same thing as in twoSum() */
vec2 twoSub(float a, float b) {
  float s = (a - b);
  float v = (s * ONE - a) * ONE;
  float err = (a - (s - v) * ONE) * ONE * ONE * ONE - (b + v);
  return vec2(s, err);
}
#else
