Re: [patch,libgfortran] Fix binary128 ERFC_SCALED

[Steve Kargl <sgk at troutmask dot apl dot washington dot edu>]

On Sun, Nov 17, 2013 at 11:29:50AM +0100, FX wrote:
> This patch fixes libgfortran?s binary128 [aka real(kind=16)] variant of ERFC_SCALED. The original code, which I had lifted from netlib, gives only 18 significant decimal digits, which is not enough for binary128 (33 decimal digits).
> 
> I thus implemented a new variant for binary128. For arguments < 12, it simply calls erfcq() then multiplies by expq(x*x). For larger arguments, it uses a power expansion in 1/x. The new implementation provides answers within to 2 ulp of the correct value.
> 
> Regtested on x86_64-apple-darwin13, comes with a testcase. OK to commit?
> FX
> 

This is probably ok as it fixes a 2**64ish ULP problem.
There is a missed optimization in

+  if (x < 12)
+    {
+      /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).
+        This is not perfect, but much better than netlib.  */
+      return erfcq(x) * expq(x*x);
+    }

If x is less than approximately -8192, then erfc(x)*exp(x*x)
overflows.  Something like the following shortcircuits the
2 function calls and their possible argument reduction steps.

  volatile GFC_REAL_16 zero = 0;   
  if (x < 12)
    {
      if (x < some_overflow_value_to_be_determined)
         return (1 / zero);
      return erfcq(x) * expq(x*x);
    }
 
PS: There is missing whitespace: x*x should be x * x.


-- 
Steve