[PATCH 23/25] math: Use erfc from CORE-MATH

[DJ Delorie]

Paul Zimmermann <Paul.Zimmermann@inria.fr> writes:
>> > +  for (int i = 0; i < 22; i++)
>> > +    if (x == exceptions[i][0])
>> > +      return exceptions[i][1] + exceptions[i][2];
>> 
>> IMHO the need for this code implies that all your other algorithms are
>> faulty, if you need to fix their errors...
>
> not necessarily. Correct rounding is ensured using Ziv's algorithm:
>
> 1) [fast path] compute an approximation y of f(x) using say 16 guard bits,
>    and an error bound e
> 2) [rounding test] if y-e and y+e round to the same value, return this value
> 3) [accurate path] using an approximation y of f(x) using a larger precision,
>    and round it to the target precision
>
> The correctness of the accurate path is validated by the knowledge of "hard-
> to-round" inputs. However, if the accurate path returns a wrong result for some
> hard-to-round input, it is sometimes more efficient to use it as an
> exceptional value, than increasing the precision of the accurate path.

Ah, ok.  The above would make a fabulous comment ;-)

>> > +  /* subnormal exceptions */
>> > +  if (x == 0x1.a8f7bfbd15495p+4)
>> > +    return fma (0x1p-1074, -0.25, 0x1.99ef5883f656cp-1024);
>> 
>> Could have been part of the exceptions table?
>
> I'm not sure, since the output should be 0x1.99ef5883f656cp-1024 for rounding
> to nearest and upward, and 0x1.99ef5883f6568p-1024 for roundings toward zero
> and downward. For the other exceptional cases, we return round(yh+yl) where
> yh+yl is a double-double approximation of erfc(x). Here since we are in the
> subnormal case, I cannot find such yh and yl values.

Perhaps its comment could read "subnormal exceptions, where the simple
math of the exception table doesn't round correctly in all cases" ?