Re: [PATCH 2/2] soft-fp: Add new KF routines

[Joseph Myers <joseph at codesourcery dot com>]

On Mon, 26 Oct 2015, Tulio Magno Quites Machado Filho wrote:

> From: Michael Meissner <meissner@linux.vnet.ibm.com>
> 
> Add cmpukf2, extendkftf2 and trunctfk2 to soft-fp.
> This is the minimal set of routines required to be imported in GCC for
> IEEE 128-bit floating point support.

As noted, these belong directly in libgcc, in an architecture-specific 
directory.  The soft-fp directory in glibc is for generic, 
architecture-independent implementations for the standard 
architecture-independent modes.  It has some files that in fact aren't 
used in glibc, but not anything inherently architecture-specific like 
this.

> +CMPtype
> +__cmpukf2 (TFtype a, TFtype b)
> +{
> +  FP_DECL_EX;
> +  FP_DECL_Q (A);
> +  FP_DECL_Q (B);
> +  CMPtype r;
> +
> +  FP_INIT_EXCEPTIONS;
> +  FP_UNPACK_RAW_Q (A, a);
> +  FP_UNPACK_RAW_Q (B, b);
> +  FP_CMP_Q (r, A, B, 2, 2);
> +  if (r == CMP_INVALID)
> +    FP_SET_EXCEPTION (FP_EX_INVALID);

No, with current FP_CMP macros you never need to set exceptions 
explicitly.  If the last argument is 2 exceptions should be raised for all 
NaN operands, 1 only for signaling NaNs (so for __cmpukf2 you want 1 as 
the last argument, for __cmpokf2 you want 2 there).

> +  FP_HANDLE_EXCEPTIONS;
> +
> +  return (r < CMP_LOW || r > CMP_HIGH) ? PPC_UNORDERED : ppc_cr_map[r-CMP_LOW];

The r value is -1, 0, 1 or 2 (the fourth argument to FP_CMP_Q).  A result 
< CMP_LOW is not possible.

> +__ibm128
> +__extendkftf2 (__float128 value)
> +{
> +  double high, low;
> +
> +  high = (double) value;
> +  if (__builtin_isnan (high) || __builtin_isinf (high))
> +    low = high;

No, that's incorrect for infinities.  See the ibm-ldouble-format file (in 
libgcc/config/rs6000/): the low part of an infinity must be 0 or -0, not 
another infinity (the low part of a NaN doesn't matter).

> +  else
> +    {
> +      low = (double) (value - (__float128)high);

There are cases where this will produce an invalid IBM long double value, 
where the low part is 0.5ulp of the high part and the high part has the 
least significant bit of the mantissa set (note that glibc contains code 
relying on this aspect of the IBM long double format (see 
sysdeps/ieee754/ldbl-128ibm/s_rintl.c: "However, if we have a canonical 
long double with the low double 0.5 or -0.5, then the high double must be 
even.")).

Thus, you need to renormalize after computing the initial approximations 
to the high and low parts.  Something like:

high_new = high + low;
low = (high - high_new) + low;
high = high_new;

> +      /* Use copysign to propigate the sign bit so that -0.0Q becomes -0.0L.  */
> +      low = __builtin_copysign (low, high);

And that's completely wrong.  The correct low part may be nonzero with 
either sign.  You need to special-case zeroes, but you can do that by 
simply including all of zero, infinity and NaN in the initial case where 
the low part is set to 0 (a low part of either +0 or -0 is fine for all 
those cases).

Also watch out for spacing in casts throughout this patch; it should be 
"(type) value" not "(type)value".

-- 
Joseph S. Myers
joseph@codesourcery.com