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
