 #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) * ONE;
-  float v = (s - a);
-  float err = (a - (s - v) * ONE) * ONE - (b + v);
+  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

  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);
  x = sub_fp64(sum_fp64(x, mul_fp64(a, exp_fp64(-x))), vec2(1.0, 0.0));
  return x;
}


  t = mul_fp64(r, INVERSE_FACTORIAL_5_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;
}
