Re: Conformance issues to address after 0.9.12 release

[Szabolcs Nagy <nsz@port70.net>]

* Rich Felker <dalias@aerifal.cx> [2013-07-29 02:34:56 -0400]:
> - i387 math asm does not truncate excess precision. Whether or not
>   this omission is conforming in terms of the return value, it results
>   in lost underflow exceptions, as demonstrated by nsz's math tests.

the underflow problem is not i387 or excess precision related:
many (odd) math functions are almost the identity function
around x==0 (sin,asin,tan,atan,sinh,atanh,..)

in those cases the traditional implementation is

  if (fabs(x) < thres)
    return x;
  if (fabs(x) < thres2)
    return x + C*x*x*x;
  ...

(in case of double precision thres is usually around 0x1p-27
so |x|*0x1p-54 > |C*x*x*x|)
(note that nan should be checked first and the thresholds are
usually checked on the bit representation with int airthmetics)

so for subnormal results x is returned: no exception is raised
(underflow and inexact should be raised if x!=0 since the
result is just an approximation)

one might think that the first check is useless optimization,
but if the C*x*x*x part is always calculated then underflow is
raised even if x is nowhere near subnormal (eg x == 0x1p-500)

a correct solution is

  if (fabs(x) < thres) {
    if (fabs(x) < 0x1p-1022)
      FORCE_EVAL(x * 1e-100); // raise inexact and underflow if x!=0
    else
      FORCE_EVAL(x + 1e100); // raise inexact
    return x;
  }
  if (fabs(x) < thres2) {
    FORCE_EVAL(x + 1e-100); // raise inexact (may be omitted)
    return x + C*x*x*x;
  }

which is ugly code bloat because volatile load/store is not
optimized in FORCE_EVAL by gcc

i387 is usually problematic in the second part: x + C*x*x*x
may not raise inexact if x is float and FLT_EVAL_METHOD==2
so the extra FORCE_EVAL or equivalent is needed