Re: [PATCHv5] New generic sinf

[Joseph Myers <joseph at codesourcery dot com>]

On Wed, 29 Nov 2017, Rajalakshmi Srinivasaraghavan wrote:

> +static const double invpio4_table[] = {
> +  0x0p+0,
> +  0x1.45f306cp+0,
> +  0x1.c9c882ap-28,
> +  0x1.4fe13a8p-58,
> +  0x1.f47d4dp-85,
> +  0x1.bb81b6cp-112,
> +  0x1.4acc9ep-142,
> +  0x1.0e4107cp-169
> +};

> +	  else                  /* |x| >= 2^23.  */
> +	    {
> +	      x = fabs (x);

Presumably this should be fabsf - no obvious reason to convert to double 
for fabs and then convert back to float.

> +	      int exponent;
> +	      GET_FLOAT_WORD (exponent, x);
> +	      exponent =
> +	        (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS;
> +	      exponent += 3;
> +	      exponent /= 28;
> +	      double a = invpio4_table[exponent] * x;
> +	      double b = invpio4_table[exponent + 1] * x;
> +	      double c = invpio4_table[exponent + 2] * x;
> +	      double d = invpio4_table[exponent + 3] * x;
> +	      unsigned long l = a;

Are you sure this is safe when unsigned long is 32-bit?  Did you test any 
32-bit configurations (you didn't specify whether your s390 testing was 
32-bit, 64-bit or both)?

Say the unbiased exponent is 52, so adding 3 and dividing by 28 gives 1.  
That is, you're computing 0x1.45f306cp+0 * x and then converting to 
unsigned long.  But with an exponent of 52, that's well outside the range 
of 32-bit long, so resulting in an INVALID exception and an unspecified 
result.  I think you need to use uint64_t here.

(Also, throughout the patch, glibc style would say "unsigned long int", 
not just "unsigned long", if you do want the machine-dependent long rather 
than needing 64-bit.)

-- 
Joseph S. Myers
joseph@codesourcery.com