math: fix asin, atan, log1p, tanh to raise underflow on subnormal

for these functions f(x)=x for small inputs, because f(0)=0 and
f'(0)=1, but for subnormal values they should raise the underflow
flag (required by annex F), if they are approximated by a polynomial
around 0 then spurious underflow should be avoided (not required by
annex F)

all these functions should raise inexact flag for small x if x!=0,
but it's not required by the standard and it does not seem a worthy
goal, so support for it is removed in some cases.

raising underflow:
- x*x may not raise underflow for subnormal x if FLT_EVAL_METHOD!=0
- x*x may raise spurious underflow for normal x if FLT_EVAL_METHOD==0
- in case of double subnormal x, store x as float
- in case of float subnormal x, store x*x as float

1 files changed, 5 insertions, 2 deletions

diff --git a/src/math/log1p.c b/src/math/log1p.c
index 6c67249c..9bed63c2 100644
--- a/src/math/log1p.c
+++ b/src/math/log1p.c
@@ -104,9 +104,12 @@ double log1p(double x)
 			return (x-x)/(x-x);         /* log1p(x<-1)=NaN */
 		}
 		if (ax < 0x3e200000) {   /* |x| < 2**-29 */
-			/* raise inexact */
-			if (two54 + x > 0.0 && ax < 0x3c900000)  /* |x| < 2**-54 */
+			/* if 0x1p-1022 <= |x| < 0x1p-54, avoid raising underflow */
+			if (ax < 0x3c900000 && ax >= 0x00100000)
 				return x;
+#if FLT_EVAL_METHOD != 0
+			FORCE_EVAL((float)x);
+#endif
 			return x - x*x*0.5;
 		}
 		if (hx > 0 || hx <= (int32_t)0xbfd2bec4) {  /* sqrt(2)/2- <= 1+x < sqrt(2)+ */