   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;
 }

  return sum_fp64(vec2(yn, 0.0), prod);
}

vec2 exp_fp64(vec2 a) {
  const float k = 512.0;
  const float inv_k = 1.0 / k;

  if (a.x <= -88.0) return vec2(0.0, 0.0);

  if (a.x == 1.0 && a.y == 0.0) return vec2(0.0, 0.0);
  if (a.x <= 0.0) return vec2(0.0 / 0.0, 0.0 / 0.0);
  vec2 x = vec2(log(a.x), 0.0);
  vec2 s;
#if defined(NVIDIA_FP64_WORKAROUND) || defined(INTEL_FP64_WORKAROUND)
  s = vec2(ONE, 0.0);
#else
  s = vec2(1.0, 0.0);
#endif

  x = sub_fp64(sum_fp64(x, mul_fp64(a, exp_fp64(-x))), s);
  return x;
}

