[PATCH v4] Fix for powerpc64 long double complex divide failure

[PATCH v4] Fix for powerpc64 long double complex divide failure

[Patrick McGehearty]

The v4 revision adds a test in libgcc/libgcc2.c for when
"__LIBGCC_TF_MANT_DIG__ == 106" to use __LIBGCC_DF_EPSILON__ instead
of __LIBGCC_TF_EPSILON__. That is specific to IBM 128-bit format long
doubles where EPSILON is very, very small and 1/EPSILON oveflows to
infinity. This change avoids the overflow without affecting any
other platform.

The revision retains the use of __LIBGCC_DF_* in .../rs6000/_divkc3.c
instead of __LIBGCC_KF_* since some older but still supported environments
do not recognize the KF floating point mode properly. This change has a
very tiny effect on the results (no cases measured in 100 million tests).

- - - -
This patch resolves the failure of powerpc64 long double complex divide
in native ibm long double format after the patch "Practical improvement
to libgcc complex divide".

The new code uses the following macros which are intended to be mapped
to appropriate values according to the underlying hardware representation.
See https://gcc.gnu.org/bugzilla/show_bug.cgi?id=101104

RBIG     a value near the maximum representation
RMIN     a value near the minimum representation
         (but not in the subnormal range)
RMIN2    a value moderately less than 1
RMINSCAL the inverse of RMIN2
RMAX2    RBIG * RMIN2  - a value to limit scaling to not overflow

When "long double" values were not using the IEEE 128-bit format but
the traditional IBM 128-bit, the previous code used the LDBL values
which caused overflow for RMINSCAL. The new code uses the DBL values.

RBIG  LDBL_MAX = 0x1.fffffffffffff800p+1022
      DBL_MAX  = 0x1.fffffffffffff000p+1022

RMIN  LDBL_MIN = 0x1.0000000000000000p-969
RMIN  DBL_MIN  = 0x1.0000000000000000p-1022

RMIN2 LDBL_EPSILON = 0x0.0000000000001000p-1022 = 0x1.0p-1074
RMIN2 DBL_EPSILON  = 0x1.0000000000000000p-52

RMINSCAL 1/LDBL_EPSILON = inf (1.0p+1074 does not fit in IBM 128-bit).
         1/DBL_EPSILON  = 0x1.0000000000000000p+52

RMAX2 = RBIG * RMIN2 = 0x1.fffffffffffff800p-52
        RBIG * RMIN2 = 0x1.fffffffffffff000p+970

The MAX and MIN values have only modest changes since the maximum and
minimum values are about the same as for double precision.  The
EPSILON field is considerably different. Due to how very small values
can be represented in the lower 64 bits of the IBM 128-bit floating
point, EPSILON is extremely small, so far beyond the desired value
that inversion of the value overflows and even without the overflow,
the RMAX2 is so small as to eliminate most usage of the test.

Instead of just replacing the use of KF_EPSILON with DF_EPSILON, we
replace all uses of KF_* with DF_*. Since the exponent fields are
essentially the same, we gain the positive benefits from the new
formula while avoiding all under/overflow issues in the #defines.

The change has been tested on gcc135.fsffrance.org and gains the
expected improvements in accuracy for long double complex divide.

libgcc/
	PR target/101104
	* config/rs6000/_divkc3.c (RBIG, RMIN, RMIN2, RMINSCAL, RMAX2):
	Use more correct values for native IBM 128-bit.
---
 libgcc/config/rs6000/_divkc3.c | 8 ++++----
 libgcc/libgcc2.c               | 5 +++++
 2 files changed, 9 insertions(+), 4 deletions(-)

diff --git a/libgcc/config/rs6000/_divkc3.c b/libgcc/config/rs6000/_divkc3.c
index a1d29d2..2b229c8 100644
--- a/libgcc/config/rs6000/_divkc3.c
+++ b/libgcc/config/rs6000/_divkc3.c
@@ -38,10 +38,10 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 #endif
 
 #ifndef __LONG_DOUBLE_IEEE128__
-#define RBIG   (__LIBGCC_KF_MAX__ / 2)
-#define RMIN   (__LIBGCC_KF_MIN__)
-#define RMIN2  (__LIBGCC_KF_EPSILON__)
-#define RMINSCAL (1 / __LIBGCC_KF_EPSILON__)
+#define RBIG   (__LIBGCC_DF_MAX__ / 2)
+#define RMIN   (__LIBGCC_DF_MIN__)
+#define RMIN2  (__LIBGCC_DF_EPSILON__)
+#define RMINSCAL (1 / __LIBGCC_DF_EPSILON__)
 #define RMAX2  (RBIG * RMIN2)
 #else
 #define RBIG   (__LIBGCC_TF_MAX__ / 2)
diff --git a/libgcc/libgcc2.c b/libgcc/libgcc2.c
index 38f935e..a0fc724 100644
--- a/libgcc/libgcc2.c
+++ b/libgcc/libgcc2.c
@@ -1906,8 +1906,13 @@ NAME (TYPE x, int m)
 # define NOTRUNC (!__LIBGCC_TF_EXCESS_PRECISION__)
 # define RBIG	(__LIBGCC_TF_MAX__ / 2)
 # define RMIN	(__LIBGCC_TF_MIN__)
+#if __LIBGCC_TF_MANT_DIG__ == 106
+# define RMIN2	(__LIBGCC_DF_EPSILON__)
+# define RMINSCAL (1 / __LIBGCC_DF_EPSILON__)
+# else
 # define RMIN2	(__LIBGCC_TF_EPSILON__)
 # define RMINSCAL (1 / __LIBGCC_TF_EPSILON__)
+#endif
 # define RMAX2	(RBIG * RMIN2)
 #else
 # error
-- 
1.8.3.1

Re: [PATCH v4] Fix for powerpc64 long double complex divide failure

[Joseph Myers]

On Thu, 26 Aug 2021, Patrick McGehearty via Gcc-patches wrote:

> The revision retains the use of __LIBGCC_DF_* in .../rs6000/_divkc3.c
> instead of __LIBGCC_KF_* since some older but still supported environments
> do not recognize the KF floating point mode properly. This change has a

That part of the patch is still wrong and should not be applied.  I don't 
know where the actual problem is in the compiler causing __LIBGCC_KF_* not 
to be defined, but using a DFmode macro there is incorrect.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: [PATCH v4] Fix for powerpc64 long double complex divide failure

[Patrick McGehearty]

I will prepare a patch without the changes to .../rs6000/_divkc3.c
but I have no way of testing it to confirm it fixes the original
complaint.


On 8/26/2021 5:09 PM, Joseph Myers wrote:
> On Thu, 26 Aug 2021, Patrick McGehearty via Gcc-patches wrote:
>
>> The revision retains the use of __LIBGCC_DF_* in .../rs6000/_divkc3.c
>> instead of __LIBGCC_KF_* since some older but still supported environments
>> do not recognize the KF floating point mode properly. This change has a
> That part of the patch is still wrong and should not be applied.  I don't
> know where the actual problem is in the compiler causing __LIBGCC_KF_* not
> to be defined, but using a DFmode macro there is incorrect.
>