[PATCH v6] Fix for powerpc64 long double complex divide failure
- - - - New in version 6: Due to an oversight (i.e. coding error), version 5 changed the use of __LIBGCC_TF_EPSILON__ to __LIBGCC_DF_EPSILON__ but not the other LIBGCC_TF values. For correct execution of the long double test case it is necessary to also switch to using __LIBGCC_DF_MIN__. For consistency we also switch to using __LIBGCC_DF_MAX__. LDBL_MIN is 2**53 times as larger than DBL_MIN. The larger value causes the code to switch the order of computation when it is not optimal, resulting in failure for one of the values in the cdivchk_ld.c test. Using DBL_MIN does not cause that failure.. There may be opportunity for further refinement of IBM128 format Long Double complex divide, but that's beyond the scope of this patch. - - - - This 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. Discussion in the patch is adjusted to reflect this limitation. It does not make any changes to .../rs6000/_divkc3.c, leaving it to use __LIBGCC_KF__*. That means the upstream gcc will not build in older IBM environments that do not recognize the KF floating point mode properly. Environments that do not need IBM longdouble support do build cleanly. - - - - This patch addresses 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) RMIN2a value moderately less than 1 RMINSCAL the inverse of RMIN2 RMAX2RBIG * 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.f800p+1022 DBL_MAX = 0x1.f000p+1022 RMIN LDBL_MIN = 0x1.p-969 RMIN DBL_MIN = 0x1.p-1022 RMIN2 LDBL_EPSILON = 0x0.1000p-1022 = 0x1.0p-1074 RMIN2 DBL_EPSILON = 0x1.p-52 RMINSCAL 1/LDBL_EPSILON = inf (1.0p+1074 does not fit in IBM 128-bit). 1/DBL_EPSILON = 0x1.p+52 RMAX2 = RBIG * RMIN2 = 0x1.f800p-52 RBIG * RMIN2 = 0x1.f000p+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. 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 * libgcc2.c (RMIN2, RMINSCAL, RMAX2): Use more correct values for native IBM 128-bit. --- libgcc/libgcc2.c | 15 +++ 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/libgcc/libgcc2.c b/libgcc/libgcc2.c index 38f935e..e66e6f0 100644 --- a/libgcc/libgcc2.c +++ b/libgcc/libgcc2.c @@ -1904,10 +1904,17 @@ NAME (TYPE x, int m) # define MODE tc # define CEXT __LIBGCC_TF_FUNC_EXT__ # define NOTRUNC (!__LIBGCC_TF_EXCESS_PRECISION__) -# define RBIG (__LIBGCC_TF_MAX__ / 2) -# define RMIN (__LIBGCC_TF_MIN__) -# define RMIN2 (__LIBGCC_TF_EPSILON__) -# define RMINSCAL (1 / __LIBGCC_TF_EPSILON__) +# if __LIBGCC_TF_MANT_DIG__ == 106 +# define RBIG (__LIBGCC_DF_MAX__ / 2) +# define RMIN (__LIBGCC_DF_MIN__) +# define RMIN2 (__LIBGCC_DF_EPSILON__) +# define RMINSCAL (1 / __LIBGCC_DF_EPSILON__) +# else +# define RBIG (__LIBGCC_TF_MAX__ / 2) +# define RMIN (__LIBGCC_TF_MIN__) +# define RMIN2(__LIBGCC_TF_EPSILON__) +# define RMINSCAL (1 / __LIBGCC_TF_EPSILON__) +# endif # define RMAX2 (RBIG * RMIN2) #else # error -- 1.8.3.1
Re: [PATCH v5] Fix for powerpc64 long double complex divide failure
Ping Note: I don't have commit privileges for gcc. On 8/27/2021 2:59 PM, Patrick McGehearty via Gcc-patches wrote: This revision (v5) 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. Discussion in the patch is adjusted to reflect this limitation. It does not make any changes to .../rs6000/_divkc3.c, leaving it to use __LIBGCC_KF__*. That means the upstream gcc will not build in older IBM environments that do not recognize the KF floating point mode properly. Environments that do not need IBM longdouble support do build cleanly. - - - - This patch addresses 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) RMIN2a value moderately less than 1 RMINSCAL the inverse of RMIN2 RMAX2RBIG * 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.f800p+1022 DBL_MAX = 0x1.f000p+1022 RMIN LDBL_MIN = 0x1.p-969 RMIN DBL_MIN = 0x1.p-1022 RMIN2 LDBL_EPSILON = 0x0.1000p-1022 = 0x1.0p-1074 RMIN2 DBL_EPSILON = 0x1.p-52 RMINSCAL 1/LDBL_EPSILON = inf (1.0p+1074 does not fit in IBM 128-bit). 1/DBL_EPSILON = 0x1.p+52 RMAX2 = RBIG * RMIN2 = 0x1.f800p-52 RBIG * RMIN2 = 0x1.f000p+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. 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 * libgcc2.c (RMIN2, RMINSCAL, RMAX2): Use more correct values for native IBM 128-bit. --- libgcc/libgcc2.c | 9 +++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/libgcc/libgcc2.c b/libgcc/libgcc2.c index 38f935e..bf45576 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__) -# define RMIN2 (__LIBGCC_TF_EPSILON__) -# define RMINSCAL (1 / __LIBGCC_TF_EPSILON__) +# 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
[PATCH v5] Fix for powerpc64 long double complex divide failure
This revision (v5) 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. Discussion in the patch is adjusted to reflect this limitation. It does not make any changes to .../rs6000/_divkc3.c, leaving it to use __LIBGCC_KF__*. That means the upstream gcc will not build in older IBM environments that do not recognize the KF floating point mode properly. Environments that do not need IBM longdouble support do build cleanly. - - - - This patch addresses 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) RMIN2a value moderately less than 1 RMINSCAL the inverse of RMIN2 RMAX2RBIG * 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.f800p+1022 DBL_MAX = 0x1.f000p+1022 RMIN LDBL_MIN = 0x1.p-969 RMIN DBL_MIN = 0x1.p-1022 RMIN2 LDBL_EPSILON = 0x0.1000p-1022 = 0x1.0p-1074 RMIN2 DBL_EPSILON = 0x1.p-52 RMINSCAL 1/LDBL_EPSILON = inf (1.0p+1074 does not fit in IBM 128-bit). 1/DBL_EPSILON = 0x1.p+52 RMAX2 = RBIG * RMIN2 = 0x1.f800p-52 RBIG * RMIN2 = 0x1.f000p+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. 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 * libgcc2.c (RMIN2, RMINSCAL, RMAX2): Use more correct values for native IBM 128-bit. --- libgcc/libgcc2.c | 9 +++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/libgcc/libgcc2.c b/libgcc/libgcc2.c index 38f935e..bf45576 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__) -# define RMIN2 (__LIBGCC_TF_EPSILON__) -# define RMINSCAL (1 / __LIBGCC_TF_EPSILON__) +# 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
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.
[PATCH v4] Fix for powerpc64 long double complex divide failure
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) RMIN2a value moderately less than 1 RMINSCAL the inverse of RMIN2 RMAX2RBIG * 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.f800p+1022 DBL_MAX = 0x1.f000p+1022 RMIN LDBL_MIN = 0x1.p-969 RMIN DBL_MIN = 0x1.p-1022 RMIN2 LDBL_EPSILON = 0x0.1000p-1022 = 0x1.0p-1074 RMIN2 DBL_EPSILON = 0x1.p-52 RMINSCAL 1/LDBL_EPSILON = inf (1.0p+1074 does not fit in IBM 128-bit). 1/DBL_EPSILON = 0x1.p+52 RMAX2 = RBIG * RMIN2 = 0x1.f800p-52 RBIG * RMIN2 = 0x1.f000p+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 v3] Fix for powerpc64 long double complex divide failure
as well as allowing the new test results to pass. I'm not sure how else to approach this issue for a short term solution. - Patrick McGehearty (patrick.mcgehea...@oracle.com) On 8/13/2021 12:12 PM, Joseph Myers wrote: On Thu, 12 Aug 2021, Patrick McGehearty via Gcc-patches wrote: If _divkc3.c is not intended to provide a version of complex divide that handles IBM-128 format, then where should that option be handled? That should be the function __divtc3. (A single libgcc binary supports multiple long double formats, so libgcc function names referring to floating-point modes need to be understood as actually referring to a particular *format*, which may or may not correspond to the named *mode* depending on the compilation options used. Thus, libgcc functions with "tf" or "tc" in their names, on configurations such as powerpc64le-linux-gnu that ever supported IBM long double, always refer to IBM long double. It's up to the back end to ensure that, when building with TFmode = binary128, TCmode and KCmode division both get mapped to __divkc3 while ICmode division gets mapped to __divtc3; when building with TFmode = IBM long double, KCmode division should still be __divkc3, ICmode division should still be __divtc3, and TCmode division should be __divtc3 in that case as well.) Do I need add a special case for #ifndef __LONG_DOUBLE_IEEE128__ in the complex divide code in libgcc/libgcc2.c? That macro is architecture-specific so shouldn't be tested there. Doing things differently if __LIBGCC_TF_MANT_DIG__ == 106 would be reasonable (if it fixes the observed problem!), however; there are a few places in generic libgcc code that check for MANT_DIG == 106 to handle IBM long double differently. And, for completeness, does gcc support LDBL for non-IEEE on any platform besides IBM? I believe TFmode is IEEE binary128 everywhere except for those powerpc configurations where it's IBM long double.
Re: [PATCH v3] Fix for powerpc64 long double complex divide failure
I see I have more to learn about gcc's interactions with IEEE-128 format vs IBM-128 format. As we discovered here, using the IBM-128 version of LDBL_EPSILON will not yield correct answers as currently coded. If _divkc3.c is not intended to provide a version of complex divide that handles IBM-128 format, then where should that option be handled? Do I need add a special case for #ifndef __LONG_DOUBLE_IEEE128__ in the complex divide code in libgcc/libgcc2.c? And, for completeness, does gcc support LDBL for non-IEEE on any platform besides IBM? - patrick On 8/12/2021 11:47 AM, Joseph Myers wrote: On Thu, 12 Aug 2021, Patrick McGehearty via Gcc-patches wrote: This file includes quad-float128.h, which does some remapping from TF to KF depending on __LONG_DOUBLE_IEEE128__. I think you probably need to have a similar __LONG_DOUBLE_IEEE128__ conditional here. If __LONG_DOUBLE_IEEE128__ is not defined, use __LIBGCC_KF_* macros instead of __LIBGCC_TF_*; if __LONG_DOUBLE_IEEE128__ is defined, use __LIBGCC_TF_* as above. (Unless the powerpc maintainers say otherwise.) - The KF version fails when in IBM128 mode while the DF version works for that mode. KFmode should always be IEEE binary128. IFmode should always be IBM long double. TFmode may be one or the other depending on command-line options. "in IBM128 mode" should mean that the compiler defaults to long double being IBM long double and TFmode being IBM long double. But in that mode, KFmode should still be IEEE binary128 and it should still be correct to use the KF constants in this file. My understanding of ibm FP mode build procedure is minimal, but it seems that the _divkc3.c routine is built for both IEEE128 and IBM128 modes. If built for IBM128 mode (i.e., compiler defaults to TFmode = IBM long double), it should still build a function __divkc3 which takes IEEE binary128 arguments and uses IEEE binary128 (KFmode) constants. If you were changing the L_divtc3 case in libgcc2.c to use different constants in the case where TFmode is IBM long double, that would make sense to me. It's changing an IEEE-only file for an IBM long double issue that doesn't make sense. If this change causes some test using IBM long double to pass where it failed before, that indicates a build system problem somewhere.
Re: [PATCH v3] Fix for powerpc64 long double complex divide failure
The key code in _divkc3.c is: #ifndef __LONG_DOUBLE_IEEE128__ #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) #define RMIN (__LIBGCC_TF_MIN__) #define RMIN2 (__LIBGCC_TF_EPSILON__) #define RMINSCAL (1 / __LIBGCC_TF_EPSILON__) #define RMAX2 (RBIG * RMIN2) #endif I added this code based on your comment of 4/20/2021: - > This file includes quad-float128.h, which does some remapping from TF to > KF depending on __LONG_DOUBLE_IEEE128__. > > I think you probably need to have a similar __LONG_DOUBLE_IEEE128__ > conditional here. If __LONG_DOUBLE_IEEE128__ is not defined, use > __LIBGCC_KF_* macros instead of __LIBGCC_TF_*; if __LONG_DOUBLE_IEEE128__ > is defined, use __LIBGCC_TF_* as above. (Unless the powerpc maintainers > say otherwise.) - The KF version fails when in IBM128 mode while the DF version works for that mode. My understanding of ibm FP mode build procedure is minimal, but it seems that the _divkc3.c routine is built for both IEEE128 and IBM128 modes. - patrick On 8/12/2021 11:19 AM, Joseph Myers wrote: On Thu, 12 Aug 2021, Patrick McGehearty via Gcc-patches wrote: 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". This description is not consistent with the patch. __divkc3 should always be using IEEE binary128 format, not IBM long double. If this code is being built for IBM long double, something is wrong somewhere else. Using the DFmode values probably makes sense for IBM long double, but not for IEEE binary128.
[PATCH v3] Fix for powerpc64 long double complex divide failure
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) RMIN2a value moderately less than 1 RMINSCAL the inverse of RMIN2 RMAX2RBIG * 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.f800p+1022 DBL_MAX = 0x1.f000p+1022 RMIN LDBL_MIN = 0x1.p-969 RMIN DBL_MIN = 0x1.p-1022 RMIN2 LDBL_EPSILON = 0x0.1000p-1022 = 0x1.0p-1074 RMIN2 DBL_EPSILON = 0x1.p-52 RMINSCAL 1/LDBL_EPSILON = inf (1.0p+1074 does not fit in IBM 128-bit). 1/DBL_EPSILON = 0x1.p+52 RMAX2 = RBIG * RMIN2 = 0x1.f800p-52 RBIG * RMIN2 = 0x1.f000p+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 1 file changed, 4 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) -- 1.8.3.1
[PATCH v2] Fix for powerpc64 long double complex divide failure
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) RMIN2a value moderately less than 1 RMINSCAL the inverse of RMIN2 RMAX2RBIG * 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.f800p+1022 DBL_MAX = 0x1.f000p+1022 RMIN LDBL_MIN = 0x1.p-969 RMIN DBL_MIN = 0x1.p-1022 RMIN2 LDBL_EPSILON = 0x0.1000p-1022 = 0x1.0p-1074 RMIN2 DBL_EPSILON = 0x1.p-52 RMINSCAL 1/LDBL_EPSILON = inf (1.0p+1074 does not fit in IBM 128-bit). 1/DBL_EPSILON = 0x1.p+52 RMAX2 = RBIG * RMIN2 = 0x1.f800p-52 RBIG * RMIN2 = 0x1.f000p+970 The MAX and MIN values have only modest changes since the exponent field for IBM 128-bit floating point values is the same size as the exponent field for IBM 64-bit floating point values. However the EPSILON field is considerably different. Due to how 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_ESPILON, 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): Fix long double complex divide for native IBM 128-bit. --- libgcc/config/rs6000/_divkc3.c | 8 1 file changed, 4 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) -- 1.8.3.1
Re: [PATCH] Fix for powerpc64 long double complex divide failure
Ping... The fix is minimal (four lines changed). I recognize that those familiar with IBM 128-bit floating point precision is a select set of people. On the plus side, tests fail without the patch and pass with the patch. - patrick On 7/8/2021 4:24 PM, Patrick McGehearty via Gcc-patches wrote: 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". See https://gcc.gnu.org/bugzilla/show_bug.cgi?id=101104 The new code uses the following macros which are intended to be mapped to appropriate values according to the underlying hardware representation. RBIG a value near the maximum representation RMIN a value near the minimum representation (but not in the subnormal range) RMIN2a value moderately less than 1 RMINSCAL the inverse of RMIN2 RMAX2RBIG * 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.f800p+1022 DBL_MAX = 0x1.f000p+1022 RMIN LDBL_MIN = 0x1.p-969 RMIN DBL_MIN = 0x1.p-1022 RMIN2 LDBL_EPSILON = 0x0.1000p-1022 = 0x1.0p-1074 RMIN2 DBL_EPSILON = 0x1.p-52 RMINSCAL 1/LDBL_EPSILON = inf (1.0p+1074 does not fit in IBM 128-bit). 1/DBL_EPSILON = 0x1.p+52 RMAX2 = RBIG * RMIN2 = 0x1.f800p-52 RBIG * RMIN2 = 0x1.f000p+970 The MAX and MIN values have only modest changes since the exponent field for IBM 128-bit floating point values is the same size as the exponent field for IBM 64-bit floating point values. However the EPSILON field is considerably different. Due to how 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. In addition, the gcc support for the KF fields (IBM native long double format) does not exist on older gcc compilers such as the default compilers on the gcc compiler farm. That adds build complexity for users who's environment is only a few years out of date. Instead of just replacing the use of KF_EPSILON with DF_ESPILON, we replace all uses of KF_* with DF_*. The change has been tested on gcc135.fsffrance.org and gains the expected improvements in accuracy for long double complex divide. libgcc/ * config/rs6000/_divkc3.c (RBIG, RMIN, RMIN2, RMINSCAL, RMAX2): Fix long double complex divide for native IBM 128-bit --- libgcc/config/rs6000/_divkc3.c | 8 1 file changed, 4 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)
[PATCH] Fix for powerpc64 long double complex divide failure
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". See https://gcc.gnu.org/bugzilla/show_bug.cgi?id=101104 The new code uses the following macros which are intended to be mapped to appropriate values according to the underlying hardware representation. RBIG a value near the maximum representation RMIN a value near the minimum representation (but not in the subnormal range) RMIN2a value moderately less than 1 RMINSCAL the inverse of RMIN2 RMAX2RBIG * 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.f800p+1022 DBL_MAX = 0x1.f000p+1022 RMIN LDBL_MIN = 0x1.p-969 RMIN DBL_MIN = 0x1.p-1022 RMIN2 LDBL_EPSILON = 0x0.1000p-1022 = 0x1.0p-1074 RMIN2 DBL_EPSILON = 0x1.p-52 RMINSCAL 1/LDBL_EPSILON = inf (1.0p+1074 does not fit in IBM 128-bit). 1/DBL_EPSILON = 0x1.p+52 RMAX2 = RBIG * RMIN2 = 0x1.f800p-52 RBIG * RMIN2 = 0x1.f000p+970 The MAX and MIN values have only modest changes since the exponent field for IBM 128-bit floating point values is the same size as the exponent field for IBM 64-bit floating point values. However the EPSILON field is considerably different. Due to how 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. In addition, the gcc support for the KF fields (IBM native long double format) does not exist on older gcc compilers such as the default compilers on the gcc compiler farm. That adds build complexity for users who's environment is only a few years out of date. Instead of just replacing the use of KF_EPSILON with DF_ESPILON, we replace all uses of KF_* with DF_*. The change has been tested on gcc135.fsffrance.org and gains the expected improvements in accuracy for long double complex divide. libgcc/ * config/rs6000/_divkc3.c (RBIG, RMIN, RMIN2, RMINSCAL, RMAX2): Fix long double complex divide for native IBM 128-bit --- libgcc/config/rs6000/_divkc3.c | 8 1 file changed, 4 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) -- 1.8.3.1
[PATCH v11] Practical improvement to libgcc complex divide
Changes in this version from Version 10: (thanks to Joseph Myers for catching these issues) In gcc/c-family/c-cppbuiltin.c Fixed three cases where XALLOCAVEC argument contained ...LIBGCC_... but should have contained ...LIBGCC__... These off by one errors managed to not fail previously due to the limited lifetime of the resulting allocation and the fact that no additional variables were placed on the stack during that lifetime. But it was an error waiting to fail at the time of some future code or compiler change. In libgcc/config/rs6000/_divkc3.c Added an #ifndef __LONG_DOUBLE_IEEE128__ and *LIBGCC_KF* values for the RBIG, RMIN, RMIN2, RMINSCAL, and RMAX2 defines. - - - - - - Correctness and performance test programs used during development of this project may be found in the attachment to: https://www.mail-archive.com/gcc-patches@gcc.gnu.org/msg254210.html Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining appropriate __LIBGCC2_* macros in c-cppbuiltin.c. Tests are added for when both parts of the denominator have exponents small enough to allow shifting any subnormal values to normal values all input values could be scaled up without risking overflow. That gained a clear improvement in accuracy. Similarly, when either numerator was subnormal and the other
Re: [PATCH v10] Practical improvement to libgcc complex divide
- ping [A sincere and special thanks for the sustained perseverance of the reviewers in pointing me in the proper direction to get this patch polished. The original proposal was June 1, 2020 and only covered double precision. The current version is dramatically better, both from extending coverage to most precisions, improving the computation for accuracy and speed, and from improving the code maintainability. - Patrick McGehearty] On 4/7/2021 3:21 PM, Patrick McGehearty via Gcc-patches wrote: Changes in this version from Version 9: Replaced all uses of alloca with XALLOCAVEC in c_cpp_builtins() in c-cppbuiltin.c Did not replace alloca elsewhere in the same file. Fixed type of name_len. Fixed prototypes for abort & exit in new test programs. Fixed spelling errors and omitted words in comments. Changed XMTYPE to AMTYPE to avoid confusion with extended types. Removed XCTYPE as unused. (A for additional) Correctness and performance test programs used during development of this project may be found in the attachment to: https://www.mail-archive.com/gcc-patches@gcc.gnu.org/msg254210.html Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining appropriate __LIBGCC2_* macros in c-cppbuiltin.c. Tests are added for when both parts of
[PATCH v10] Practical improvement to libgcc complex divide
Changes in this version from Version 9: Replaced all uses of alloca with XALLOCAVEC in c_cpp_builtins() in c-cppbuiltin.c Did not replace alloca elsewhere in the same file. Fixed type of name_len. Fixed prototypes for abort & exit in new test programs. Fixed spelling errors and omitted words in comments. Changed XMTYPE to AMTYPE to avoid confusion with extended types. Removed XCTYPE as unused. (A for additional) Correctness and performance test programs used during development of this project may be found in the attachment to: https://www.mail-archive.com/gcc-patches@gcc.gnu.org/msg254210.html Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining appropriate __LIBGCC2_* macros in c-cppbuiltin.c. Tests are added for when both parts of the denominator have exponents small enough to allow shifting any subnormal values to normal values all input values could be scaled up without risking overflow. That gained a clear improvement in accuracy. Similarly, when either numerator was subnormal and the other numerator and both denominator values were not too large, scaling could be used to reduce risk of computing with subnormals. The test and scaling values used all fit within the allowed exponent range for each precision required by the C standard. Float precision has more
Re: [PATCH v9] Practical improvement to libgcc complex divide
I would very much like to get all the details right on my next submission.
I'm submitting this query for feedback before I do the full submission.
Most changes requested require obvious fixes. One has me confused about
the best action to take.
- - - -
On the issue of alloca vs XALLOCAVEC, I discovered that XALLOCAVEC has
been around and recommended since at least 2010, but I could not find
any discussion of why XALLOCAVEC is preferred to alloca. Understanding
the purpose of making the change would help in resolve the following
decision.
Currently,in gcc/c-family/c-cppbuiltin.c, XALLOCAVEC does not appear
at all. alloca appears 13 times. Both appear in other files in the
same directory. In the routine I'm modifying, four current allocations
to macro_name use alloca.
I see four options:
1) Use alloca as that is consistent with the existing style within
the source file. That is my current choice, but disliked by the
reviewer.
2) If I use XALLOCAVEC just for my three new uses of alloca, which
also assign to macro_name, we will have an inconsistent coding
pattern. That violates the principle of maintaining existing style
within a source file or routine. This choice seems ugly.
3) If I use XALLOCATE for my three new uses for alloca and
change the other four allcoations to macro_name in the same
routine but leave the other 9 uses of alloca alone, we have
consistency in the modified routine but not the source file.
If this is the recommended procedure, I'm willing.
4) If I change all uses of alloc to XALLOCATE in the file, that makes
the file consistent but adds source file changes to the patch in
routines that are completely unrelated to my goal of improving complex
divide. It increases risk of causing subtle problems either through
my editing errors or lack of understanding by me of any differences
between alloca and XALLOCATE. I'm less comfortable with this option.
- - - -
- - - -
The rest appear to be editing issues on my part.
I'll will make the following changes.
- - - -
I will fix the type of name_len to const size_t.
- - - -
I will fix abort/exit in the three test files to be:
extern void abort(void);
extern void exit(int status);
- - - -
I will proofread and fix the comments in the three test files.
than/that, missing "to", etc.
- - - -
In libgcc/config/rs6000/_divkc3.c
I see one additional blank line that will be removed.
- - - -
In libgcc/libgcc2.c
- - - -
I admit the names XCTYPE and XMTYPE are not very imaginative. My
intention was represent "extra precision" but the potential confusion
with extended precision did not occur to me at the time. The
XCTYPE/XMTYPE is intended to be the next larger precision to allow the
same source code to be used for both the half precision and float
precision cases.
Maybe it would be better to think "additional precision",
and use ACTYPE and AMTYPE?
- - - -
Comment "...errors than can..." should be "...errors that can..."
Will make the change.
On 3/31/2021 12:03 PM, Bernhard Reutner-Fischer wrote:
On Fri, 26 Mar 2021 23:14:41 +
Patrick McGehearty via Gcc-patches wrote:
diff --git a/gcc/c-family/c-cppbuiltin.c b/gcc/c-family/c-cppbuiltin.c
index 9f993c4..c0d9e57 100644
--- a/gcc/c-family/c-cppbuiltin.c
+++ b/gcc/c-family/c-cppbuiltin.c
@@ -1277,8 +1277,10 @@ c_cpp_builtins (cpp_reader *pfile)
{
scalar_float_mode mode = mode_iter.require ();
const char *name = GET_MODE_NAME (mode);
+ const int name_len = strlen (name);
strlen returns a size_t
Funny that we do not have a pre-seeded mode_name_len array but call
strlen on modes over and over again.
+ char float_h_prefix[16] = "";
char *macro_name
- = (char *) alloca (strlen (name)
+ = (char *) alloca (name_len
+ sizeof ("__LIBGCC__MANT_DIG__"));
sprintf (macro_name, "__LIBGCC_%s_MANT_DIG__", name);
builtin_define_with_int_value (macro_name,
@@ -1286,20 +1288,29 @@ c_cpp_builtins (cpp_reader *pfile)
if (!targetm.scalar_mode_supported_p (mode)
|| !targetm.libgcc_floating_mode_supported_p (mode))
continue;
- macro_name = (char *) alloca (strlen (name)
+ macro_name = (char *) alloca (name_len
+ sizeof ("__LIBGCC_HAS__MODE__"));
sprintf (macro_name, "__LIBGCC_HAS_%s_MODE__", name);
cpp_define (pfile, macro_name);
- macro_name = (char *) alloca (strlen (name)
+ macro_name = (char *) alloca (name_len
+ sizeof ("__LIBGCC__FUNC_EXT__"));
sprintf (macro_name, "__LIBGCC_%s_FUNC_EXT__", name);
The above should have been split out as separate independent patchlet
to get it out of the way.
As noted already the use of alloca is a pre-existing (coding-style) bug
and we should
Re: [PATCH v9] Practical improvement to libgcc complex divide
Thank you for your interest in this project.
On 3/27/2021 6:18 PM, Bernhard Reutner-Fischer wrote:
On Fri, 26 Mar 2021 23:14:41 +
Patrick McGehearty via Gcc-patches wrote:
Changes in Version 9 since Version 8:
Revised code to meet gcc coding standards in all files, especially
with respect to adding spaces around operations and removing
excess () in #define macro definitions.
Did you gather cycle counter diffs for current against new for the
normal and some edge cases?
I don't see additional value in looking at cycle counters in addition
to the timing information I presented. I prefer wall clock time to
cycle counters, perhaps because cycle counters are less reliably
implemented on some platforms (that's a separate discussion, though).
Every platform will have different counters and different timing.
Let me run through a rough operation count analysis of the double
precision case. I will also postulate instruction latencies
that might be typical of a recent processor. Real processors
will different somewhat but not by large amounts.
The current method uses Smith's method. Showing half of that method:
if (fabs(c) < fabs(d))
{
ratio = c / d;
denom = (c * ratio) + d;
xx = ((a * ratio) + b) / denom;
yy = ((b * ratio) - a) / denom;
}
It requires:
two fabs computations (typically 1 cycle each)
one floating point compare (typically 2 cycles)
1 branch decision (50-50 prediction reliabilty)
(1 cycle if predicted, 9 cycles if not predicted; ave=5)
3 multiplys, 3 divides, 3 add/subtract operations
(assume 5 cycle latency for add/mul and 15 cycle for divide)
(75 cycles)
Total: 84 cycles (typical average, old method)
I omit the NAN analysis as that is unchanged.
The cost of the new method is data dependent.
First, assume all scaling tests are false and ratio does not underflow.
Further, assume that data is consistently in those ranges, to allow
all branch predictions except the first to be predicted.
I'll also assume we have an excellent branch predictor mechanism
that can handle predicting
Then we have:
5 additional fabs, and 4 additional compare and branches (all predicted)
and one || operations for 18 additional cycles.
[This optimistic count assumes fabs(a,b) > RMIN, allowing
us to bypass the fabs < RMAX2 tests.
New total: 84+18 = 102 (optimistic average, new method)
If we can't bypass the fabs < RMAX2 tests, we need to add
four more fabs, and four && operations, and four compare
operations for 12 more cycles.
Possible total: 84+18+12 = 114
Finally, if we hit a scaling operation, then we have four
multiply operations for an additional 20 cycles plus
we've mispredicted the branch for 9 more.
114+28 = 144 cycles.
This analysis is not particularly accurate. It does not account for
some architectures having several floating point units or for
speculative execution hiding latency or many other variations in
computer architecture. Still, its good enough for a ballpark estimate
that says the typical case will be around 15-25% slower and the most
extreme cases will be around 70% slower.
When we look at the patch writeup, I report that for double precision,
the "moderate set" which represents the typical case, we see 4% to 24%
measured slowdown. The "full set" which represents having a
substantial number of values which might cause branch predictors to go
wrong as costing 38% to 56% overhead. Apparently my 'back of the
envelope' analysis above is not too far off.
Oh, I failed to cover the case where fabs(ratio) < RMIN.
That happens in the full set about 1-2% of the time.
It requires two additional divide operations for an extra
30 cycles, plus maybe 10 cycles for the mispredicted branch.
That might bump up the worst case to be about double the
normal case. That's acceptable when you consider that slow
path prevents the code from returning a completely wrong result.
As a quick aside, the test for ratio underflowing reduces the error
rate by a factor of 4 from 1.70% to 0.35% in the full exponent range
data set with about half the total overhead of the new method. Adding
all the other tests and scaling code reduces errors larger than 12
bits to roughly 1 in 10 million and 18 bits or larger to roughly 1 in
100 million. Those remaining errors are due to the issue of
subtraction cancellation, which is a harder problem than avoiding
underflow/overflow issues.
Can you please share data for -Os between current and your proposed new?
As for data for -Os for current and proposed new, it's not very interesting.
There are no differences as complex divide is not inlined by gcc.
Instead, there is a simple call to __divdc3.
I checked for inlining with gcc 4.8, gcc8.3, and gcc11.0.1 as well as with
a compiler that has my code changes.
If one uses -fcx-fortran-rules then Smith's method for complex divide is
inlined with -Os.
If one uses -fcx-limited-range then the naive method for complex divide
is inli
[PATCH v9] Practical improvement to libgcc complex divide
Changes in Version 9 since Version 8: Revised code to meet gcc coding standards in all files, especially with respect to adding spaces around operations and removing excess () in #define macro definitions. Major revision to gcc/testsuite/gcc.c-torture/execute/ieee/cdivchkld.c Prior code was focused on x86 80-bit implementation of long doubles. The new version provides tests for: IEEE 128-bit format, 80-bit format, ibm extended precision format (11-bit exponent, 113 bit mantissa), and when long double is treated as IEEE 64-bit doubles (also 11-bit exponent). The limits tested are now based on LDBL_MIN and LDBL_MANT_DIG which automatically adapts to the appropriate format. The IEEE 128-bit and 80-bit formats share input test values due to having matching 15-bit exponent sizes as doe the IBM extended double and normal 64-bit double having 11-bit exponent sizes. The input values are accurate to the lesser mantissa with zero fill for the extended length mantissas. When the test with the smaller mantissa is active, the expected result will automatically be truncated to match the available mantissa. The size of the exponent (LDBL_MAX_EXP) is used to determine which values are tested. The program was tested using x86 (80-bit), aarch64 (128-bit), IBM extended format and IEEE 64-bit format. Some adjustments where made due to IBM extended format having slightly less range at the bottom end than IEEE 64-bit format in spite of having the same size exponent field. For all four cases, the cdivchkld.c program will abort() under the old complex divide and exit() under the new complex divide. Correctness and performance test programs used during development of this project may be found in the attachment to: https://www.mail-archive.com/gcc-patches@gcc.gnu.org/msg254210.html Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor
Re: [PATCH v8] Practical improvement to libgcc complex divide
Ping... On 2/22/2021 5:08 PM, Patrick McGehearty via Gcc-patches wrote: Changes in this version from Version 7: gcc/c-family/c-cppbuiltin.c Changed modename to float_h_prefix Removed (mode == TYPE_MODE...) code left over from earlier approach. libgcc/libgcc2.c Removed conditional use of XF constants in TF case. Also left over from earlier approach. Tested _Float64x case on x86. Works as expected. Correctness and performance test programs used during development of this project may be found in the attachment to: https://www.mail-archive.com/gcc-patches@gcc.gnu.org/msg254210.html Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining appropriate __LIBGCC2_* macros in c-cppbuiltin.c. Tests are added for when both parts of the denominator have exponents small enough to allow shifting any subnormal values to normal values all input values could be scaled up without risking overflow. That gained a clear improvement in accuracy. Similarly, when either numerator was subnormal and the other numerator and both denominator values were not too large, scaling could be used to reduce risk of computing with subnormals. The test and scaling values used all fit within the allowed exponent range for e
[PATCH v8] Practical improvement to libgcc complex divide
Changes in this version from Version 7: gcc/c-family/c-cppbuiltin.c Changed modename to float_h_prefix Removed (mode == TYPE_MODE...) code left over from earlier approach. libgcc/libgcc2.c Removed conditional use of XF constants in TF case. Also left over from earlier approach. Tested _Float64x case on x86. Works as expected. Correctness and performance test programs used during development of this project may be found in the attachment to: https://www.mail-archive.com/gcc-patches@gcc.gnu.org/msg254210.html Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining appropriate __LIBGCC2_* macros in c-cppbuiltin.c. Tests are added for when both parts of the denominator have exponents small enough to allow shifting any subnormal values to normal values all input values could be scaled up without risking overflow. That gained a clear improvement in accuracy. Similarly, when either numerator was subnormal and the other numerator and both denominator values were not too large, scaling could be used to reduce risk of computing with subnormals. The test and scaling values used all fit within the allowed exponent range for each precision required by the C standard. Float precision has more difficulty with getting correct answers than double
Re: [PATCH v7] Practical improvement to libgcc complex divide
I decided to use float_h_prefix instead of modename or mode prefix as Joseph has better insight as to naming expectations when others are working in gcc/c-family/c-cppbuiltin.c. All requested changes made and tested. Patch coming shortly. - patrick On 2/19/2021 5:48 PM, Joseph Myers wrote: On Fri, 19 Feb 2021, Patrick McGehearty via Gcc-patches wrote: Here you're properly computing the mapping from mode to float.h macro prefix (though I think "modename" is a confusing name for the variable used for float.h macro prefixes; to me, "modename" suggests the names such as SF or DF, which are actually "name" here, and something like "float_h_prefix" would be better than "modename"). float_h_prefix puts me in mind of half-precision floating point rather than float.h. How would modeprefix do instead of modename? I'm not sure modeprefix is better; you could say SF is the prefix of SFmode, for example. Hence my suggestion of a name that makes clear it's the prefix *in float.h*. typeprefix might be a possibility, as it's really a prefix associated with a corresponding type rather than directly with the mode.
Re: [PATCH v7] Practical improvement to libgcc complex divide
Comments inline
On 2/19/2021 3:27 PM, Joseph Myers wrote:
On Tue, 2 Feb 2021, Patrick McGehearty via Gcc-patches wrote:
if (mode == TYPE_MODE (double_type_node))
- ; /* Empty suffix correct. */
+ {
+ ; /* Empty suffix correct. */
+ memcpy (modename, "DBL", 4);
+ }
else if (mode == TYPE_MODE (float_type_node))
- suffix[0] = 'f';
+ {
+ suffix[0] = 'f';
+ memcpy (modename, "FLT", 4);
+ }
else if (mode == TYPE_MODE (long_double_type_node))
- suffix[0] = 'l';
+ {
+ suffix[0] = 'l';
+ memcpy (modename, "LDBL", 4);
+ }
else
{
bool found_suffix = false;
@@ -1305,6 +1316,8 @@ c_cpp_builtins (cpp_reader *pfile)
sprintf (suffix, "f%d%s", floatn_nx_types[i].n,
floatn_nx_types[i].extended ? "x" : "");
found_suffix = true;
+ sprintf (modename, "FLT%d%s", floatn_nx_types[i].n,
+floatn_nx_types[i].extended ? "X" : "");
break;
}
gcc_assert (found_suffix);
Here you're properly computing the mapping from mode to float.h macro
prefix (though I think "modename" is a confusing name for the variable
used for float.h macro prefixes; to me, "modename" suggests the names such
as SF or DF, which are actually "name" here, and something like
"float_h_prefix" would be better than "modename").
float_h_prefix puts me in mind of half-precision floating point rather
than float.h.
How would modeprefix do instead of modename?
+ if ((mode == TYPE_MODE (float_type_node))
+ || (mode == TYPE_MODE (double_type_node))
+ || (mode == TYPE_MODE (long_double_type_node)))
But then here you still have the check for the mode matching one of those
three types, which defeats the point of computing the prefix in a way that
works for *all* floating-point modes supported in libgcc. (The above
"gcc_assert (found_suffix)" asserts that the relevant mapping did get
found.)
To be able to use the __LIBGCC_* macros in libgcc in all cases, they
should be defined, using the modename (or float_h_prefix) found above,
whether or not the mode matches one of float, double and long double.
My apologies. You are completely correct. I intended to remove the "mode
== TYPE_MODE" code
when making the other changes but overlooked it.
#elif defined(L_multc3) || defined(L_divtc3)
# define MTYPETFtype
# define CTYPETCtype
# define MODE tc
# define CEXT __LIBGCC_TF_FUNC_EXT__
# define NOTRUNC (!__LIBGCC_TF_EXCESS_PRECISION__)
+#if defined(__LIBGCC_TF_MIN__)
+# define RBIG ((__LIBGCC_TF_MAX__)/2)
+# define RMIN (__LIBGCC_TF_MIN__)
+# define RMIN2 (__LIBGCC_TF_EPSILON__)
+# define RMINSCAL (1/__LIBGCC_TF_EPSILON__)
+#else
+# define RBIG ((__LIBGCC_XF_MAX__)/2)
+# define RMIN (__LIBGCC_XF_MIN__)
+# define RMIN2 (__LIBGCC_XF_EPSILON__)
+# define RMINSCAL (1/__LIBGCC_XF_EPSILON__)
+#endif
Once you've removed the unnecessary conditional above, you don't need this
conditional on __LIBGCC_TF_MIN__ being defined either; the *TF* macros
will always be defined when TFmode is supported in libgcc, so you never
need to use *XF* macros instead when building TFmode code.
And as you note, the test for LIBGCC_TF_MIN being defined is not needed
either.
I'll resubmit the patch after I rebuild and complete testings on the
revised code.
- patrick
Re: [PATCH v7] Practical improvement to libgcc complex divide
Ping - Submitted Feb 2, 2021 I believe this version fixes all issues raised in previous submissions. - Patrick McGehearty On 2/2/2021 4:25 PM, Patrick McGehearty via Gcc-patches wrote: Changes in this version from Version 6: Updated copyrights for following three files to -2021. gcc/c-family/c-cppbuiltin.c Moved code for setting LIBGCC modename to FUNC_EXT section. Added code to set modename for any additional modes such as FLT128 or FLT64X. libgcc/libgcc2.c Changed RMIN2 & RMINSCAL for XF and TF modes to use the matching LIBGCC_?F_EPSILON instead of LIBGCC_DF_EPSILON libgcc/config/rs6000/_divkc3.c Changed RMIN2 & RMINSCAL for TF mode to use the matching LIBGCC_TF_EPSILON instead of LIBGCC_DF_EPSILON Resubmitted with correction to following link. No other changes from previous PATCH v7. Correctness and performance test programs used during development of this project may be found in the attachment to: https://gcc.gnu.org/pipermail/gcc-patches/2021-January/564239.html with the attachment at https://gcc.gnu.org/pipermail/gcc-patches/attachments/20210125/b5a7df3b/attachment-0001.bin Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining
Re: [PATCH v6] Practical improvement to libgcc complex divide
I saw the gcc.gnu.org at the end of the address and did not realize the link was not part gnu.org. The gcc.gnu.org link is: https://gcc.gnu.org/pipermail/gcc-patches/2021-January/564239.html with the attachment at https://gcc.gnu.org/pipermail/gcc-patches/attachments/20210125/b5a7df3b/attachment-0001.bin I'll resubmit the patch with the corrected link. - patrick On 2/1/2021 5:49 PM, Joseph Myers wrote: On Mon, 1 Feb 2021, Patrick McGehearty via Gcc-patches wrote: The message which contains the attached gzipped tarball of the development test programs is: https://www.mail-archive.com/gcc-patches@gcc.gnu.org/msg254210.html I'll include that link in the new patch as well. I think it's best to include the URL in the archives on gcc.gnu.org rather than on a third-party site.
[PATCH v7] Practical improvement to libgcc complex divide
Changes in this version from Version 6: Updated copyrights for following three files to -2021. gcc/c-family/c-cppbuiltin.c Moved code for setting LIBGCC modename to FUNC_EXT section. Added code to set modename for any additional modes such as FLT128 or FLT64X. libgcc/libgcc2.c Changed RMIN2 & RMINSCAL for XF and TF modes to use the matching LIBGCC_?F_EPSILON instead of LIBGCC_DF_EPSILON libgcc/config/rs6000/_divkc3.c Changed RMIN2 & RMINSCAL for TF mode to use the matching LIBGCC_TF_EPSILON instead of LIBGCC_DF_EPSILON Resubmitted with correction to following link. No other changes from previous PATCH v7. Correctness and performance test programs used during development of this project may be found in the attachment to: https://gcc.gnu.org/pipermail/gcc-patches/2021-January/564239.html with the attachment at https://gcc.gnu.org/pipermail/gcc-patches/attachments/20210125/b5a7df3b/attachment-0001.bin Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining appropriate __LIBGCC2_* macros in c-cppbuiltin.c. Tests are added for when both parts of the denominator have exponents small enough to allow shifting any subnormal values to normal values all input values could be scaled up without risking
[PATCH v7] Practical improvement to libgcc complex divide
Changes in this version from Version 6: Updated copyrights for following three files to -2021. gcc/c-family/c-cppbuiltin.c Moved code for setting LIBGCC modename to FUNC_EXT section. Added code to set modename for any additional modes such as FLT128 or FLT64X. libgcc/libgcc2.c Changed RMIN2 & RMINSCAL for XF and TF modes to use the matching LIBGCC_?F_EPSILON instead of LIBGCC_DF_EPSILON libgcc/config/rs6000/_divkc3.c Changed RMIN2 & RMINSCAL for TF mode to use the matching LIBGCC_TF_EPSILON instead of LIBGCC_DF_EPSILON Correctness and performance test programs used during development of this project may be found in the attachment to: https://www.mail-archive.com/gcc-patches@gcc.gnu.org/msg254210.html Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining appropriate __LIBGCC2_* macros in c-cppbuiltin.c. Tests are added for when both parts of the denominator have exponents small enough to allow shifting any subnormal values to normal values all input values could be scaled up without risking overflow. That gained a clear improvement in accuracy. Similarly, when either numerator was subnormal and the other numerator and both denominator values were not too large, scaling could be used to reduce
Re: [PATCH v6] Practical improvement to libgcc complex divide
The message which contains the attached gzipped tarball of the
development test programs is:
https://www.mail-archive.com/gcc-patches@gcc.gnu.org/msg254210.html
I'll include that link in the new patch as well.
I can't recall why I did not use _LIBGCC_XF_EPSILON__ or
_LIBGCC_TF_EPSILON__ before. Perhaps a copy/paste oversight?
Fixed now.
On the issue of handling more float_modes than the well known FLT,
DBL, and LDBL, it took me some time to understand how and where the
additional modes were defined. For example, it is hard to find FLT64
when you grep for FLOAT64. Once I found gcc/glinclude/float.h, and
some other files, everything fell into place. Actually making the code
changes as suggested to support all float modes was easy.
The new patch [v7] will be forwarded as soon as testing is complete,
likely later today.
- patrick
On 1/21/2021 5:04 PM, Joseph Myers wrote:
On Fri, 18 Dec 2020, Patrick McGehearty via Gcc-patches wrote:
TEST Data
I'd still like to see the test data / code used to produce the accuracy
and performance results made available somewhere (presumably with a link
then being provided in the commit message).
+ if ((mode == TYPE_MODE (float_type_node))
+ || (mode == TYPE_MODE (double_type_node))
+ || (mode == TYPE_MODE (long_double_type_node)))
+ {
+ char val_name[64];
+ char fname[8] = "";
+ if (mode == TYPE_MODE (float_type_node))
+ memcpy (fname, "FLT", 4);
+ else if (mode == TYPE_MODE (double_type_node))
+ memcpy (fname, "DBL", 4);
+ else if (mode == TYPE_MODE (long_double_type_node))
+ memcpy (fname, "LDBL", 5);
This logic being used to generate EPSILON, MAX and MIN macros only handles
modes that match float, double or long double (so won't define the macros
for a mode that only matches another type such as _Float128, for example).
Earlier in the same function, there is existing code to define
__LIBGCC_%s_FUNC_EXT__. That code has to do something similar, to
determine the matching type for a mode - but it also handles _FloatN /
_FloatNx types, and has an assertion at the end that some matching type
was found.
Rather than having this code which handles a more limited set of types, I
think the __LIBGCC_%s_FUNC_EXT__ code should be extended, so that as well
as computing a function suffix it also computes a prefix such as FLT, DBL,
FLT128 or FLT64X. Then all supported floating-point modes can get these
three LIBGCC macros defined, rather than just those matching float, double
or long double.
#elif defined(L_mulxc3) || defined(L_divxc3)
# define MTYPEXFtype
# define CTYPEXCtype
# define MODE xc
# define CEXT __LIBGCC_XF_FUNC_EXT__
# define NOTRUNC (!__LIBGCC_XF_EXCESS_PRECISION__)
+# define RBIG ((__LIBGCC_XF_MAX__)/2)
+# define RMIN (__LIBGCC_XF_MIN__)
+# define RMIN2 (__LIBGCC_DF_EPSILON__)
+# define RMINSCAL (1/__LIBGCC_DF_EPSILON__)
+# define RMAX2 ((RBIG)*(RMIN2))
I'd then expect __LIBGCC_XF_EPSILON__ to be used for XFmode in place of
__LIBGCC_DF_EPSILON__ unless there is some good reason to use the latter
(which would need a comment to explain it if so).
#elif defined(L_multc3) || defined(L_divtc3)
# define MTYPETFtype
# define CTYPETCtype
# define MODE tc
# define CEXT __LIBGCC_TF_FUNC_EXT__
# define NOTRUNC (!__LIBGCC_TF_EXCESS_PRECISION__)
+#if defined(__LIBGCC_TF_MIN__)
+# define RBIG ((__LIBGCC_TF_MAX__)/2)
+# define RMIN (__LIBGCC_TF_MIN__)
+#else
+# define RBIG ((__LIBGCC_XF_MAX__)/2)
+# define RMIN (__LIBGCC_XF_MIN__)
+#endif
+# define RMIN2 (__LIBGCC_DF_EPSILON__)
+# define RMINSCAL (1/__LIBGCC_DF_EPSILON__)
And, likewise, with the suggested changes to c-cppbuiltin.c this code can
use __LIBGCC_TF_MAX__, __LIBGCC_TF_MIN__ and __LIBGCC_TF_EPSILON__
unconditionally, without ever needing to use XF or DF macros. (If you
want to use a different EPSILON value in the case where TFmode is IBM long
double because of LDBL_EPSILON being too small in that case, condition
that on __LIBGCC_TF_MANT_DIG__ == 106, and use ((TFtype) 0x1p-106) in
place of __LIBGCC_TF_EPSILON__ in that case.)
Re: [PATCH v6] Practical improvement to libgcc complex divide
There are two open gcc bugzilla bugs for complex divide. 59714 complex divide is surprising on targets with FMA 60646 investigate improved complex division algorithms Either might be appropriate. 59714 started my investigation, and my solution resolves that issue, while my solution seems closer to the request of 60646, except 60646 seems to target the middle-end instead of libgcc. Anyway, once I created the tarball and gzipped it, I found it shrank to 100Kbytes. Rather than have further delays, I've attached the gzip file to this email. If someone has trouble unzipping, I can email the .tar file to them directly (500Kbytes). - patrick mcgehearty On 1/25/2021 2:48 PM, Joseph Myers wrote: On Mon, 25 Jan 2021, Patrick McGehearty wrote: On 1/21/2021 5:04 PM, Joseph Myers wrote: On Fri, 18 Dec 2020, Patrick McGehearty via Gcc-patches wrote: TEST Data I'd still like to see the test data / code used to produce the accuracy and performance results made available somewhere (presumably with a link then being provided in the commit message). I've packaged up the code/scripts/README and some data in a directory. I had it ready last time and meant to provide it, but it slipped my mind while I was preparing the patch. Total size is under a megabyte, so moving it around should not be a problem. Is this sort of development support code usually kept somewhere in the gcc tree or should it be placed elsewhere, like in my personal github directory? I think that sort of thing (when it's too big to just attach to the email with the patch - if it can be attached to the email, or to a relevant open bug in Bugzilla, that's the simplest thing to do) would typically go somewhere else rather than in the GCC tree. cdiv_measure.tar.gz Description: GNU Zip compressed data
Re: [PATCH v6] Practical improvement to libgcc complex divide
On 1/21/2021 5:04 PM, Joseph Myers wrote:
On Fri, 18 Dec 2020, Patrick McGehearty via Gcc-patches wrote:
TEST Data
I'd still like to see the test data / code used to produce the accuracy
and performance results made available somewhere (presumably with a link
then being provided in the commit message).
I've packaged up the code/scripts/README and some data in a directory.
I had it ready last time and meant to provide it, but it slipped my mind
while I was preparing the patch. Total size is under a megabyte, so moving
it around should not be a problem.
Is this sort of development support code usually kept somewhere in the
gcc tree or
should it be placed elsewhere, like in my personal github directory?
- patrick
+ if ((mode == TYPE_MODE (float_type_node))
+ || (mode == TYPE_MODE (double_type_node))
+ || (mode == TYPE_MODE (long_double_type_node)))
+ {
+ char val_name[64];
+ char fname[8] = "";
+ if (mode == TYPE_MODE (float_type_node))
+ memcpy (fname, "FLT", 4);
+ else if (mode == TYPE_MODE (double_type_node))
+ memcpy (fname, "DBL", 4);
+ else if (mode == TYPE_MODE (long_double_type_node))
+ memcpy (fname, "LDBL", 5);
This logic being used to generate EPSILON, MAX and MIN macros only handles
modes that match float, double or long double (so won't define the macros
for a mode that only matches another type such as _Float128, for example).
Earlier in the same function, there is existing code to define
__LIBGCC_%s_FUNC_EXT__. That code has to do something similar, to
determine the matching type for a mode - but it also handles _FloatN /
_FloatNx types, and has an assertion at the end that some matching type
was found.
Rather than having this code which handles a more limited set of types, I
think the __LIBGCC_%s_FUNC_EXT__ code should be extended, so that as well
as computing a function suffix it also computes a prefix such as FLT, DBL,
FLT128 or FLT64X. Then all supported floating-point modes can get these
three LIBGCC macros defined, rather than just those matching float, double
or long double.
#elif defined(L_mulxc3) || defined(L_divxc3)
# define MTYPEXFtype
# define CTYPEXCtype
# define MODE xc
# define CEXT __LIBGCC_XF_FUNC_EXT__
# define NOTRUNC (!__LIBGCC_XF_EXCESS_PRECISION__)
+# define RBIG ((__LIBGCC_XF_MAX__)/2)
+# define RMIN (__LIBGCC_XF_MIN__)
+# define RMIN2 (__LIBGCC_DF_EPSILON__)
+# define RMINSCAL (1/__LIBGCC_DF_EPSILON__)
+# define RMAX2 ((RBIG)*(RMIN2))
I'd then expect __LIBGCC_XF_EPSILON__ to be used for XFmode in place of
__LIBGCC_DF_EPSILON__ unless there is some good reason to use the latter
(which would need a comment to explain it if so).
#elif defined(L_multc3) || defined(L_divtc3)
# define MTYPETFtype
# define CTYPETCtype
# define MODE tc
# define CEXT __LIBGCC_TF_FUNC_EXT__
# define NOTRUNC (!__LIBGCC_TF_EXCESS_PRECISION__)
+#if defined(__LIBGCC_TF_MIN__)
+# define RBIG ((__LIBGCC_TF_MAX__)/2)
+# define RMIN (__LIBGCC_TF_MIN__)
+#else
+# define RBIG ((__LIBGCC_XF_MAX__)/2)
+# define RMIN (__LIBGCC_XF_MIN__)
+#endif
+# define RMIN2 (__LIBGCC_DF_EPSILON__)
+# define RMINSCAL (1/__LIBGCC_DF_EPSILON__)
And, likewise, with the suggested changes to c-cppbuiltin.c this code can
use __LIBGCC_TF_MAX__, __LIBGCC_TF_MIN__ and __LIBGCC_TF_EPSILON__
unconditionally, without ever needing to use XF or DF macros. (If you
want to use a different EPSILON value in the case where TFmode is IBM long
double because of LDBL_EPSILON being too small in that case, condition
that on __LIBGCC_TF_MANT_DIG__ == 106, and use ((TFtype) 0x1p-106) in
place of __LIBGCC_TF_EPSILON__ in that case.)
Re: [PATCH v6] Practical improvement to libgcc complex divide
Ping. With the holidays coming shortly after I posted this this on Dec 17, 2020, I understand any delays in getting it reviewed. Still, I thought I should give a ping. I very much would like to complete this project while I can still find all my notes related to the design decisions. :-) - Patrick McGehearty (patrick.mcgehea...@oracle.com) On 12/17/2020 9:43 PM, Patrick McGehearty via Gcc-patches wrote: Changes in this version from Version 5: Commit message change: Information for Changelogs put in proper format. Thanks go to Jakub Jelinek for enlightening me how to do that and providing an excellent example. c-cppbuiltin.c - Fixed various formatting issues regarding misplaced or missing spaces. Replaced strncpy with memcpy. Replaced all uses of strlen (name) with name_len in calls to alloca. Did not replace alloca with XALLOCAVEC. There are no other uses of XALLOCAVEC in this file, although XALLOCAVEC was used in c-ada-spec.c and c-common.c. Following the principle of least change consistent coding practice within a single source file I am hesitant to make the change. If gcc maintainers desire to generally replace all uses of alloca with XALLOCVEC at this time for this file, I'm willing to make the change. No changes to the other modified files from Version 5. Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for
[PATCH v6] Practical improvement to libgcc complex divide
Changes in this version from Version 5: Commit message change: Information for Changelogs put in proper format. Thanks go to Jakub Jelinek for enlightening me how to do that and providing an excellent example. c-cppbuiltin.c - Fixed various formatting issues regarding misplaced or missing spaces. Replaced strncpy with memcpy. Replaced all uses of strlen (name) with name_len in calls to alloca. Did not replace alloca with XALLOCAVEC. There are no other uses of XALLOCAVEC in this file, although XALLOCAVEC was used in c-ada-spec.c and c-common.c. Following the principle of least change consistent coding practice within a single source file I am hesitant to make the change. If gcc maintainers desire to generally replace all uses of alloca with XALLOCVEC at this time for this file, I'm willing to make the change. No changes to the other modified files from Version 5. Summary of Purpose This patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test data by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining appropriate __LIBGCC2_* macros in c-cppbuiltin.c. Tests are added for when both parts of the denominator have exponents small enough to allow shifting any subnormal values to normal values all input values could be scaled up without risking
Re: [PATCH v5] Practical Improvement to libgcc Complex Divide
On 12/10/2020 5:35 PM, Jakub Jelinek wrote: On Thu, Dec 10, 2020 at 10:27:46AM -0600, Patrick McGehearty via Gcc-patches wrote: Thank you for your rapid feedback. I'll fix the various formatting issues (spaces in the wrong places and such as well as revise the Changelog magic) in the next submission. It will wait for Joseph's review to also make any changes he suggests. I'll also try to train myself to be more sensitive to gcc formatting conventions while proofreading. I'm reluctant to change or use XALLOCAVEC instead of alloca as that is not the current style elsewhere in the routine. If so, I can fix it incrementally. But, at least fix up the formatting, that was the reason I've mentioned XALLOCAVEC, because the alloca call formatting was off. If there is a stated policy somewhere that says "replace alloca with XALLOCA* whenever a file is to be updated/changed, I'll follow that, but so far I have not found that statement. Without the statement, iIf it's a good idea to change alloca to XALLOCAVEC, I'd support that as a separate patch. Currently there are 25 uses of alloca in gcc/c-family/* and 6 uses of XALLOCA*. There are more other uses of alloca in the gcc src tree if you are looking to make a clean sweep of alloca. I was not reading this list when XALLOCA* first started getting used (quick googling suggests it was a decade ago) so I'm not up on why it makes a difference. No need to confuse the issues by doing some changes with the complex divide patch and some changes separately. As we don't know when complex divide will make it through review, if you get XALLOCAVEC changes committed to c-cppbuiltin.c before complex divide is approved, I will have no problem with replacing alloca with XALLOCAVEC in my patch. And I will certainly do my best to fix all formatting issues in either case. On the strcpy, strncpy, and memcpy question, given short length of the string being copied, I don't think it makes much difference. The two other copy operations in the file are memcpy. memcpy might be slightly better since it is generally more frequently seen and more likely that gcc has special case code to inline short fixed length memcpy as a few assignments. Even if both strncpy and memcpy are inlined, the memcpy code may be simplier as it does not need to be concerned with special treatment of nulls. I'll change the strncpy to memcpy. Even if short, strncpy is a badly designed API that in 99% cases just shouldn't be used. Either the string is shorter than the passed argument, then in most cases completely useless zeroes the rest of the buffer (exception is security sensitive code that needs to overwrite everything that has been before), or it is the same length and then one should use strcpy instead, or there is string truncation which doesn't zero terminate, which is very rarely something one wants to do. I generally avoid strcpy for various reasons that may not apply in this particular case. I used strncpy simply because "I need to copy a string, what does that? oh, strncpy." I was not sensitive to the points you make about strncpy. memcpy is clearly right for this use case. Jakub - patrick
Re: [PATCH v5] Practical Improvement to libgcc Complex Divide
Thank you for your rapid feedback.
I'll fix the various formatting issues (spaces in the wrong places
and such as well as revise the Changelog magic) in the next submission.
It will wait for Joseph's review to also make any changes he suggests.
I'll also try to train myself to be more sensitive to gcc formatting
conventions while proofreading.
I'm reluctant to change or use XALLOCAVEC instead of alloca as that
is not the current style elsewhere in the routine.
I agree that namelen=strlen(name) would be an apparent optimization,
but since *name is declared const char, I would think the compiler
would need to compute strlen(name) only one time. Again, I'm reluctant
to change the existing code patterns.
On the strcpy, strncpy, and memcpy question, given short length of
the string being copied, I don't think it makes much difference.
The two other copy operations in the file are memcpy.
memcpy might be slightly better since it is generally more frequently
seen and more likely that gcc has special case code to inline
short fixed length memcpy as a few assignments. Even if both strncpy
and memcpy are inlined, the memcpy code may be simplier as it does
not need to be concerned with special treatment of nulls.
I'll change the strncpy to memcpy.
- patrick
On 12/8/2020 5:31 PM, Jakub Jelinek wrote:
On Tue, Dec 08, 2020 at 10:32:33PM +, Patrick McGehearty via Gcc-patches
wrote:
2020-12-08 Patrick McGehearty
* gcc/c-family/c-cppbuiltin.c - Add supporting macros for new complex divide.
* libgcc2.c (__divsc3, __divdc3, __divxc3, __divtc3): Improve complex divide.
* libgcc/config/rs6000/_divkc3.c - Complex divide changes for rs6000.
* gcc/testsuite/gcc.c-torture/execute/ieee/cdivchkd.c - double cdiv test.
* gcc/testsuite/gcc.c-torture/execute/ieee/cdivchkf.c - float cdiv test.
* gcc/testsuite/gcc.c-torture/execute/ieee/cdivchkld.c - long double cdiv test.
Thanks for working on this, I'll defer review to Joseph, just want to add a few
random comments.
The above ChangeLog will not get through the commit checking scripts,
one needs two spaces before and after name instead of just one,
pathnames should be relative to the corresponding ChangeLog file and one
should separate what goes to each ChangeLog, and lines except the first one
should be tab indented. So it should look like:
2020-12-08 Patrick McGehearty
gcc/c-family/
* c-cppbuiltin.c (c_cpp_builtins): Add supporting macros for new
complex divide.
libgcc/
* libgcc2.c (XMTYPE, XCTYPE, RBIG, RMIN, RMIN2, RMINSCAL, RMAX2):
Define.
(__divsc3, __divdc3, __divxc3, __divtc3): Improve complex divide.
* config/rs6000/_divkc3.c (RBIG, RMIN, RMIN2, RMINSCAL, RMAX2):
Define.
(__divkc3): Improve complex divide.
gcc/testsuite/
* gcc.c-torture/execute/ieee/cdivchkd.c: New test.
* gcc.c-torture/execute/ieee/cdivchkf.c: New test.
* gcc.c-torture/execute/ieee/cdivchkld.c: New test.
or so.
--- a/gcc/c-family/c-cppbuiltin.c
+++ b/gcc/c-family/c-cppbuiltin.c
@@ -1347,6 +1347,47 @@ c_cpp_builtins (cpp_reader *pfile)
"PRECISION__"));
sprintf (macro_name, "__LIBGCC_%s_EXCESS_PRECISION__", name);
builtin_define_with_int_value (macro_name, excess_precision);
+
+ if ((mode == TYPE_MODE (float_type_node))
+ || (mode == TYPE_MODE (double_type_node))
+ || (mode == TYPE_MODE (long_double_type_node)))
+ {
+ char val_name[64];
+ char fname[8] = "";
+ if (mode == TYPE_MODE (float_type_node))
+ strncpy(fname, "FLT",4);
Formatting, there should be space before ( for calls, and space in between
, and 4. Also, what is the point of using strncpy? strcpy or
memcpy would do.
+ else if (mode == TYPE_MODE (double_type_node))
+ strncpy(fname, "DBL",4);
+ else if (mode == TYPE_MODE (long_double_type_node))
+ strncpy(fname, "LDBL",5);
+
+ if ( (mode == TYPE_MODE (float_type_node))
+ || (mode == TYPE_MODE (double_type_node)) )
Formatting, no spaces in between the ( ( and ) ).
+ {
+ macro_name = (char *) alloca (strlen (name)
+ + sizeof ("__LIBGCC_EPSILON__"
+ ));
This should use XALLOCAVEC macro, so
macro_name
= XALLOCAVEC (char, strlen (name)
+ sizeof ("__LIBGCC_EPSILON__"));
I admit it is a preexisting problem in the code above it too.
+ sprintf (macro_name, "__LIBGCC_%s_EPSILON__", name);
+ sprintf( val_name, "__%s_EPSILON__", fname);
Space before ( rather than after it.
+ builtin_define_with_val
[PATCH v5] Practical Improvement to libgcc Complex Divide
Summary of Purpose The following patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test set by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rare and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Further investigation showed changing the half precision algorithm to use the simple formula (real=a*c+b*d imag=b*c-a*d) caused no loss of precision and modest improvement in performance. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. This use is made generic by defining appropriate __LIBGCC2_* macros in c-cppbuiltin.c. Tests are added for when both parts of the denominator have exponents small enough to allow shifting any subnormal values to normal values all input values could be scaled up without risking overflow. That gained a clear improvement in accuracy. Similarly, when either numerator was subnormal and the other numerator and both denominator values were not too large, scaling could be used to reduce risk of computing with subnormals. The test and scaling values used all fit within the allowed exponent range for each precision required by the C standard. Float precision has more difficulty with getting correct answers than double precision. When hardware for double precision floating point operations is available, float precision is now handled in double precision intermediate calculations with the simple algorithm the same as the half-precision method of using float precision for intermediate calculations. Using the higher precision yields exact results for all tested input values (64-bit double, 32-bit float) with the only performance cost being the requirement to convert the four input values from float to double. If double precision hardware is not available, then
Re: [PATCH] Practical Improvement to libgcc Complex Divide
It took some work, but I think I've responded to all the issues raised here.
Patch V5 coming right after this mail.
On 11/16/2020 8:34 PM, Joseph Myers wrote:
On Tue, 8 Sep 2020, Patrick McGehearty via Gcc-patches wrote:
This project started with an investigation related to
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1]
provided an overview of past and recent practice for computing complex
divide. The current glibc implementation is based on Robert Smith's
algorithm [2] from 1962. A google search found the paper by Baudin
and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's
proposed patch [4] is based on that paper.
Thanks, I've now read Baudin and Smith so can review the patch properly.
I'm fine with the overall algorithm, so my comments generally relate to
how the code should best be integrated into libgcc while keeping it
properly machine-mode-generic as far as possible.
I developed two sets of test set by randomly distributing values over
a restricted range and the full range of input values. The current
Are these tests available somewhere?
After some polishing, the development tests are now ready to share.
I've got them in a single directory (a README, 47 mostly small .c files,
various scripts for running tests and sample outputs from all the tests.
the tarball totals about 0.5 MBytes.) These tests are intended for
developing and comparing different complex divide algorithms.
They are NOT intended or structured for routine compiler testing.
The complex divide compiler tests are in the accompanying patch
and are discussed later in this note.
Support for half, float, double, extended, and long double precision
is included as all are handled with suitable preprocessor symbols in a
single source routine. Since half precision is computed with float
precision as per current libgcc practice, the enhanced algorithm
provides no benefit for half precision and would cost performance.
Therefore half precision is left unchanged.
The existing constants for each precision:
float: FLT_MAX, FLT_MIN;
double: DBL_MAX, DBL_MIN;
extended and/or long double: LDBL_MAX, LDBL_MIN
are used for avoiding the more common overflow/underflow cases.
In general, libgcc code works with modes, not types; hardcoding references
to a particular mapping between modes and types is problematic. Rather,
the existing code in c-cppbuiltin.c that has a loop over modes should be
extended to provide whatever information is needed, as macros defined for
each machine mode.
/* For libgcc-internal use only. */
if (flag_building_libgcc)
{
/* Properties of floating-point modes for libgcc2.c. */
opt_scalar_float_mode mode_iter;
FOR_EACH_MODE_IN_CLASS (mode_iter, MODE_FLOAT)
{
[...]
For example, that defines macros such as __LIBGCC_DF_FUNC_EXT__ and
__LIBGCC_DF_MANT_DIG__. The _FUNC_EXT__ definition involves that code
computing a mapping to types.
I'd suggest defining additional macros such as __LIBGCC_DF_MAX__ in the
same code - for each supported floating-point mode. They can be defined
to __FLT_MAX__ etc. (the predefined macros rather than the ones in
float.h) - the existing code that computes a suffix for functions can be
adjusted so it also computes the string such as "FLT", "DBL", "LDBL",
"FLT128" etc.
(I suggest defining to __FLT_MAX__ rather than to the expansion of
__FLT_MAX__ because that avoids any tricky interactions with the logic to
compute such expansions lazily. I suggest __FLT_MAX__ rather than the
FLT_MAX name from float.h because that way you avoid any need to define
feature test macros to access names such as FLT128_MAX.)
After some study, I've done my best to follow your recommendations
for using modes. I've defined __LIBGCC_xx_yyy__, where xx is SF, DF,
XF, TF and yyy is MIN, MAX, and EPSILON in c-cppbuiltin.c. SF uses FLT,
DF used DBL, and XF and TF use LDBL. There is no need for those values
in HF mode because the HF code always uses SF precision.
diff --git a/gcc/c-family/c-cppbuiltin.c b/gcc/c-family/c-cppbuiltin.c
index 74ecca8..02c06d8 100644
--- a/gcc/c-family/c-cppbuiltin.c
+++ b/gcc/c-family/c-cppbuiltin.c
@@ -1343,6 +1343,11 @@ c_cpp_builtins (cpp_reader *pfile)
builtin_define_with_value ("__LIBGCC_INIT_SECTION_ASM_OP__",
INIT_SECTION_ASM_OP, 1);
#endif
+ /* For libgcc float/double optimization */
+#ifdef HAVE_adddf3
+ builtin_define_with_int_value ("__LIBGCC_HAVE_HWDBL__",
+HAVE_adddf3);
+#endif
This is another thing to handle more generically - possibly with something
like the mode_has_fma function, and defining a macro for each mode, named
after the mode, rather than only for DFmode. For an alternative, see the
discussion below.
Defining this value generically but not using/testing it seems
more likely to be subject it future issues when someone tries
to
Re: [PATCH] Practical Improvement to libgcc Complex Divide
Joseph, thank you for your detailed review and comments.
I will get to work on the necessary revisions as well
as find for a suitable place for sharing my random number
generating tests.
- patrick
On 11/16/2020 8:34 PM, Joseph Myers wrote:
On Tue, 8 Sep 2020, Patrick McGehearty via Gcc-patches wrote:
This project started with an investigation related to
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1]
provided an overview of past and recent practice for computing complex
divide. The current glibc implementation is based on Robert Smith's
algorithm [2] from 1962. A google search found the paper by Baudin
and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's
proposed patch [4] is based on that paper.
Thanks, I've now read Baudin and Smith so can review the patch properly.
I'm fine with the overall algorithm, so my comments generally relate to
how the code should best be integrated into libgcc while keeping it
properly machine-mode-generic as far as possible.
I developed two sets of test set by randomly distributing values over
a restricted range and the full range of input values. The current
Are these tests available somewhere?
Support for half, float, double, extended, and long double precision
is included as all are handled with suitable preprocessor symbols in a
single source routine. Since half precision is computed with float
precision as per current libgcc practice, the enhanced algorithm
provides no benefit for half precision and would cost performance.
Therefore half precision is left unchanged.
The existing constants for each precision:
float: FLT_MAX, FLT_MIN;
double: DBL_MAX, DBL_MIN;
extended and/or long double: LDBL_MAX, LDBL_MIN
are used for avoiding the more common overflow/underflow cases.
In general, libgcc code works with modes, not types; hardcoding references
to a particular mapping between modes and types is problematic. Rather,
the existing code in c-cppbuiltin.c that has a loop over modes should be
extended to provide whatever information is needed, as macros defined for
each machine mode.
/* For libgcc-internal use only. */
if (flag_building_libgcc)
{
/* Properties of floating-point modes for libgcc2.c. */
opt_scalar_float_mode mode_iter;
FOR_EACH_MODE_IN_CLASS (mode_iter, MODE_FLOAT)
{
[...]
For example, that defines macros such as __LIBGCC_DF_FUNC_EXT__ and
__LIBGCC_DF_MANT_DIG__. The _FUNC_EXT__ definition involves that code
computing a mapping to types.
I'd suggest defining additional macros such as __LIBGCC_DF_MAX__ in the
same code - for each supported floating-point mode. They can be defined
to __FLT_MAX__ etc. (the predefined macros rather than the ones in
float.h) - the existing code that computes a suffix for functions can be
adjusted so it also computes the string such as "FLT", "DBL", "LDBL",
"FLT128" etc.
(I suggest defining to __FLT_MAX__ rather than to the expansion of
__FLT_MAX__ because that avoids any tricky interactions with the logic to
compute such expansions lazily. I suggest __FLT_MAX__ rather than the
FLT_MAX name from float.h because that way you avoid any need to define
feature test macros to access names such as FLT128_MAX.)
diff --git a/gcc/c-family/c-cppbuiltin.c b/gcc/c-family/c-cppbuiltin.c
index 74ecca8..02c06d8 100644
--- a/gcc/c-family/c-cppbuiltin.c
+++ b/gcc/c-family/c-cppbuiltin.c
@@ -1343,6 +1343,11 @@ c_cpp_builtins (cpp_reader *pfile)
builtin_define_with_value ("__LIBGCC_INIT_SECTION_ASM_OP__",
INIT_SECTION_ASM_OP, 1);
#endif
+ /* For libgcc float/double optimization */
+#ifdef HAVE_adddf3
+ builtin_define_with_int_value ("__LIBGCC_HAVE_HWDBL__",
+HAVE_adddf3);
+#endif
This is another thing to handle more generically - possibly with something
like the mode_has_fma function, and defining a macro for each mode, named
after the mode, rather than only for DFmode. For an alternative, see the
discussion below.
diff --git a/libgcc/ChangeLog b/libgcc/ChangeLog
index ccfd6f6..8bd66c5 100644
--- a/libgcc/ChangeLog
+++ b/libgcc/ChangeLog
@@ -1,3 +1,10 @@
+2020-08-27 Patrick McGehearty
+
+ * libgcc2.c (__divsc3, __divdc3, __divxc3, __divtc3): Enhance
+ accuracy of complex divide by avoiding underflow/overflow when
+ ratio underflows or when arguments have very large or very
+ small exponents.
Note that diffs to ChangeLog files should not now be included in patches;
the ChangeLog content needs to be included in the proposed commit message
instead for automatic ChangeLog generation.
+#if defined(L_divsc3)
+#define RBIG ((FLT_MAX)/2.0)
+#define RMIN (FLT_MIN)
+#define RMIN2 (0x1.0p-21)
+#define RMINSCAL (0x1.0p+19)
+#define RMAX2 ((RBIG)*(RMIN2))
+#endif
I'd expect all of these to use generic macros based on the mode. For the
division by 2.0, probably also divide by in
Re: [PATCH] Practical Improvement to libgcc Complex Divide
2nd ping - submitted Sept 8, 2020, last comment Sept 9, 2020 This patch is a correctness bug fix and not a change to any user visible interface. It makes a major improvement to complex divide accuracy, including algorithm improvements known for 8 years. It is highly desirable to get these changes into gcc11 and not wait another year. In recent years there have been substantial improvements in mathematical functions in glibc and gcc. I'm hoping to continue those contributions with the goal of making glibc/gcc the "go to" solution for portable mathematical software development. Substantial delays in getting reviews makes it more difficult to get management backing for additional mathematical improvement projects. Also, if reviewer has any unexpected questions about design decisions, it would be helpful to ask them while I still remember why I made those choices. - Patrick McGehearty patrick.mcgehea...@oracle.com On 10/13/2020 9:54 AM, Patrick McGehearty via Gcc-patches wrote: Ping - still need review of version 4 of this patch. It has been over a month since the last comment. - patrick On 9/9/2020 2:13 AM, Richard Biener wrote: On Tue, Sep 8, 2020 at 8:50 PM Patrick McGehearty via Gcc-patches wrote: (Version 4) (Added in version 4) Fixed Changelog entry to include __divsc3, __divdc3, __divxc3, __divtc3. Revised description to avoid incorrect use of "ulp (units last place)". Modified float precison case to use double precision when double precision hardware is available. Otherwise float uses the new algorithm. Added code to scale subnormal numerator arguments when appropriate. This change reduces 16 bit errors in double precision by a factor of 140. Revised results charts to match current version of code. Added background of tuning approach. Summary of Purpose The following patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://urldefense.com/v3/__https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714__;!!GqivPVa7Brio!NjjEhnQQ38VNyP_v8nlAm9uVjvZldUnobfY5hZdq22cMMVauop64MFw3nOHIQXUmy8PToRw$ . Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test set by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rarest and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, do
Re: [PATCH] Practical Improvement to libgcc Complex Divide
Ping - still need review of version 4 of this patch. It has been over a month since the last comment. - patrick On 9/9/2020 2:13 AM, Richard Biener wrote: On Tue, Sep 8, 2020 at 8:50 PM Patrick McGehearty via Gcc-patches wrote: (Version 4) (Added in version 4) Fixed Changelog entry to include __divsc3, __divdc3, __divxc3, __divtc3. Revised description to avoid incorrect use of "ulp (units last place)". Modified float precison case to use double precision when double precision hardware is available. Otherwise float uses the new algorithm. Added code to scale subnormal numerator arguments when appropriate. This change reduces 16 bit errors in double precision by a factor of 140. Revised results charts to match current version of code. Added background of tuning approach. Summary of Purpose The following patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://urldefense.com/v3/__https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714__;!!GqivPVa7Brio!NjjEhnQQ38VNyP_v8nlAm9uVjvZldUnobfY5hZdq22cMMVauop64MFw3nOHIQXUmy8PToRw$ . Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test set by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rarest and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Therefore half precision is left unchanged. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. Testing for when both parts of the denominator had exponents roughly small enough to allow shifting any subnormal values to normal values, all input values could be scaled up without risking unnecessary overflow and gaining a clear improvement in accuracy. Similarly, when either numerator was subnormal and the other numerator and both denominator values were not too large, scaling could be used to reduce risk of computing with subnormals. The test
Re: [PATCH] Practical Improvement to libgcc Complex Divide
On 9/9/2020 2:13 AM, Richard Biener wrote: Thanks for working on this. Speaking about performance and accuracy I spot a few opportunities to use FMAs [and eventually vectorization] - do FMAs change anything on the accuracy analysis (is there the chance they'd make it worse?). We might want to use IFUNCs in libgcc to version for ISA variants (with/without FMA)? Thanks, Richard. Richard, Thank you for bringing up the issue of fused multiply-add (fma). All the results I presented in the latest patch were measured with fma active. That's because in my early testing I ran experiments and found that fma was consistently more accurate than no fma. In response to your query, I repeated that set of tests on my final submission and present them in the following table. Number of results out of 10 million with greater than or equal to the listed number of bits in error. full range limited exponents no fma with fma no fma with fma 1 bits= 20088 16664 34479 24707 2 bits= 1110 900 2359 1762 3 bits= 518 440 1163 882 4 bits= 197 143 612 445 5 bits= 102 72 313 232 6 bits= 49 43 170 119 7 bits= 25 21 82 49 8 bits= 16 11 33 26 9 bits= 9 5 14 14 10 bits= 3 3 8 4 11 bits= 2 2 3 2 12 bits= 1 1 0 2 No differences for 13 or greater bits. Errors for both cases drop off rapidly as we increase the number of bits required for a result to be considered an error. While using fma shows a consistent advantage in fewer errors, there are cases were no fma gives a more accurate answer. A detailed examination of the full range/7 bit case which is listed as having 25 errors greater than 7 bits for "no fma" and 21 errors greater than 7 bits for "fma" In that test, 1 case had the same size error for both 8 cases had a larger error with no fma 21 cases had a larger error with fma. Further examination showed those differences were generally less than two bits. That summary makes clear that while using fma does not always give a better answer, it occasionally provides a slight improvement. Even in the limited exponent case where 1 bit difference is counted as an error, using fma or not using fma only shows a different result about 1 time in 1000. While interesting, that size and frequency of difference is not enough to support having two versions the library routine in my mind. I'd rather put further effort into improving the accuracy or performance of other libgcc/glibc math functions. Just as one example, Paul Zimmerman's work shows opportunities. For more details, see: https://urldefense.com/v3/__https://members.loria.fr/PZimmermann/papers/accuracy.pdf__;!!GqivPVa7Brio!PREWxi54-6JnIBbz8jjKEYGoZ3x6Nz5_4dXoalIf8uR1i3NKHHCgdGZJbzEXQmRMrmKmk38$ - patrick
[PATCH] Practical Improvement to libgcc Complex Divide
(Version 4) (Added in version 4) Fixed Changelog entry to include __divsc3, __divdc3, __divxc3, __divtc3. Revised description to avoid incorrect use of "ulp (units last place)". Modified float precison case to use double precision when double precision hardware is available. Otherwise float uses the new algorithm. Added code to scale subnormal numerator arguments when appropriate. This change reduces 16 bit errors in double precision by a factor of 140. Revised results charts to match current version of code. Added background of tuning approach. Summary of Purpose The following patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [2] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. Background This project started with an investigation related to https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714. Study of Beebe[1] provided an overview of past and recent practice for computing complex divide. The current glibc implementation is based on Robert Smith's algorithm [2] from 1962. A google search found the paper by Baudin and Smith [3] (same Robert Smith) published in 2012. Elen Kalda's proposed patch [4] is based on that paper. I developed two sets of test set by randomly distributing values over a restricted range and the full range of input values. The current complex divide handled the restricted range well enough, but failed on the full range more than 1% of the time. Baudin and Smith's primary test for "ratio" equals zero reduced the cases with 16 or more error bits by a factor of 5, but still left too many flawed answers. Adding debug print out to cases with substantial errors allowed me to see the intermediate calculations for test values that failed. I noted that for many of the failures, "ratio" was a subnormal. Changing the "ratio" test from check for zero to check for subnormal reduced the 16 bit error rate by another factor of 12. This single modified test provides the greatest benefit for the least cost, but the percentage of cases with greater than 16 bit errors (double precision data) is still greater than 0.027% (2.7 in 10,000). Continued examination of remaining errors and their intermediate computations led to the various tests of input value tests and scaling to avoid under/overflow. The current patch does not handle some of the rarest and most extreme combinations of input values, but the random test data is only showing 1 case in 10 million that has an error of greater than 12 bits. That case has 18 bits of error and is due to subtraction cancellation. These results are significantly better than the results reported by Baudin and Smith. Support for half, float, double, extended, and long double precision is included as all are handled with suitable preprocessor symbols in a single source routine. Since half precision is computed with float precision as per current libgcc practice, the enhanced algorithm provides no benefit for half precision and would cost performance. Therefore half precision is left unchanged. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. Testing for when both parts of the denominator had exponents roughly small enough to allow shifting any subnormal values to normal values, all input values could be scaled up without risking unnecessary overflow and gaining a clear improvement in accuracy. Similarly, when either numerator was subnormal and the other numerator and both denominator values were not too large, scaling could be used to reduce risk of computing with subnormals. The test and scaling values used all fit within the allowed exponent range for each precision required by the C standard. Float precision has even more difficulty with getting correct answers than double precision. When hardware for double precision floating point operations is available, float precision is now handled in double precision intermediate calculations with the original Smith algorithm (i.e. the
Re: PING [Patch][Middle-end]Add -fzero-call-used-regs=[skip|used-gpr|all-gpr|used|all]
My understanding is this feature/flag is not intended to be "default on". It is intended to be used in security sensitive environments such as the Linux kernel where it was requested by kernel security experts. I'm not understanding the objection here if the feature is requested by security teams and the average cost is modest. My background is in performance and application optimization. I agree that for typical computation oriented, non-secure applications, I would not use the feature, but for system applications that have the ability to cross protection boundaries, it seems to be clearly a worthwhile feature. - patrick On 9/7/2020 9:44 AM, Segher Boessenkool wrote: On Fri, Sep 04, 2020 at 01:23:14AM +, Rodriguez Bahena, Victor wrote: Qing, thanks a lot for the measurement, I am not sure if this is the limit of overhead the community is willing to accept by adding extra security (me as gcc user will be willing to accept). The overhead is of course bearable for most programs / users, but what is the return? For what percentage of programs are ROP attacks no longer possible, for example. Segher
Re: [PATCH V3] Practical Improvement to libgcc Complex Divide
2nd ping. Any estimate on when a reviewer might get to this improvement in the accuracy of Complex Divide? I'm happy to supply more info on what testing I've done and details about design decisions. I'd prefer to do that sooner than later as who knows when corporate priority decisions might prevent me from having time for rapid response. - Patrick McGehearty (patrick.mcgehea...@oracle.com) On 7/21/2020 12:19 PM, Patrick McGehearty via Gcc-patches wrote: Ping On 7/1/2020 11:30 AM, Patrick McGehearty via Gcc-patches wrote: (Version 3) (Added in version 3) Support for half, float, extended, and long double precision has been added to the prior work for double precision. Since half precision is computed with float precision as per current libgcc practice, the enhanced underflow/overflow tests provide no benefit for half precision and would cost performance. Therefore half precision is left unchanged. The existing constants for each precision: float: FLT_MAX, FLT_MIN; double: DBL_MAX, DBL_MIN; extended and/or long double: LDBL_MAX, LDBL_MIN are used for avoiding the more common overflow/underflow cases. Additional investigation showed that testing for when both parts of the denominator had exponents roughly small enough to allow shifting any subnormal values to normal values, all input values could be scaled up without risking unnecessary overflow and gaining a clear improvement in accuracy. The test and scaling values used all fit within the allowed exponent range for each precision required by the C standard. The remaining number of troubling results in version 3 is measurably smaller than in versions 1 and 2. The timing and precision tables below have been revised appropriately to match the algorithms used in this version for double precision and additional tables added to include results for other precisions. In prior versions, I omitted mention of the bug report that started me on this project: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714 complex division is surprising on targets with FMA (was: on aarch64) With the proposed method, whether using FMA or not, dividing 1.0+3.0i by 1.0+3.0i correctly returns 1.0+0.0i. I also have added a reference to Beebe's "The Mathematical Function Computation Handbook" [4] which was my starting point for research into better complex divide methods. (Added for Version 2) In my initial research, I missed Elen Kalda's proposed patch https://gcc.gnu.org/legacy-ml/gcc-patches/2019-08/msg01629.html [3] Thanks to Joseph Myers for providing me with the pointer. This version includes performance and accuracy comparisions between Elen's proposed patch and my latest patch version for double precision. (from earlier Versions) The following patch to libgcc/libgcc2.c __divdc3 provides an opportunity to gain important improvements to the quality of answers for the default complex divide routine (half, float, double, extended, long double precisions) when dealing with very large or very small exponents. The current code correctly implements Smith's method (1962) [1] further modified by c99's requirements for dealing with NaN (not a number) results. When working with input values where the exponents are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are substantially different from the answers provided by quad precision more than 1% of the time. This error rate may be unacceptable for many applications that cannot a priori restrict their computations to the safe range. The proposed method reduces the frequency of "substantially different" answers by more than 99% for double precision at a modest cost of performance. Differences between current gcc methods and the new method will be described. Then accuracy and performance differences will be discussed. NOTATION For all of the following, the notation is: Input complex values: a+bi (a= real part, b= imaginary part) c+di Output complex value: e+fi = (a+bi)/(c+di) For the result tables: current = current method (SMITH) b1div = method proposed by Elen Kalda b2div = alternate method considered by Elen Kalda new1 = new method using 1 divide and 2 multiplies new = new method proposed by this patch DESCRIPTIONS of different complex divide methods: NAIVE COMPUTATION (-fcx-limited-range): e = (a*c + b*d)/(c*c + d*d) f = (b*c - a*d)/(c*c + d*d) Note that c*c and d*d will overflow or underflow if either c or d is outside the range 2^-538 to 2^512. This method is available in gcc when the switch -fcx-limited-range is used. That switch is also enabled by -ffast-math. Only one who has a clear understanding of the maximum range of intermediate values generated by a computation should consider using this switch. SMITH's METHOD (current libgcc): if(fabs(c) /* Changing the order of operations avoids the underflow of r impacting the result. */ x = (a + (d * (b / c))) * t; y = (b - (d * (a / c))) * t; } } if
Re: [PATCH V3] Practical Improvement to libgcc Complex Divide
Ping
On 7/1/2020 11:30 AM, Patrick McGehearty via Gcc-patches wrote:
(Version 3)
(Added in version 3)
Support for half, float, extended, and long double precision has
been added to the prior work for double precision. Since half precision
is computed with float precision as per current libgcc practice,
the enhanced underflow/overflow tests provide no benefit for half
precision and would cost performance. Therefore half precision is
left unchanged.
The existing constants for each precision:
float: FLT_MAX, FLT_MIN;
double: DBL_MAX, DBL_MIN;
extended and/or long double: LDBL_MAX, LDBL_MIN
are used for avoiding the more common overflow/underflow cases.
Additional investigation showed that testing for when both parts of
the denominator had exponents roughly small enough to allow shifting
any subnormal values to normal values, all input values could be
scaled up without risking unnecessary overflow and gaining a clear
improvement in accuracy. The test and scaling values used all fit
within the allowed exponent range for each precision required by the C
standard. The remaining number of troubling results in version 3 is
measurably smaller than in versions 1 and 2.
The timing and precision tables below have been revised appropriately
to match the algorithms used in this version for double precision
and additional tables added to include results for other precisions.
In prior versions, I omitted mention of the bug report that started me
on this project: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714
complex division is surprising on targets with FMA (was: on aarch64)
With the proposed method, whether using FMA or not, dividing
1.0+3.0i by 1.0+3.0i correctly returns 1.0+0.0i.
I also have added a reference to Beebe's "The Mathematical Function
Computation Handbook" [4] which was my starting point for research
into better complex divide methods.
(Added for Version 2)
In my initial research, I missed Elen Kalda's proposed patch
https://gcc.gnu.org/legacy-ml/gcc-patches/2019-08/msg01629.html [3]
Thanks to Joseph Myers for providing me with the pointer.
This version includes performance and accuracy comparisions
between Elen's proposed patch and my latest patch version for
double precision.
(from earlier Versions)
The following patch to libgcc/libgcc2.c __divdc3 provides an
opportunity to gain important improvements to the quality of answers
for the default complex divide routine (half, float, double, extended,
long double precisions) when dealing with very large or very small exponents.
The current code correctly implements Smith's method (1962) [1]
further modified by c99's requirements for dealing with NaN (not a
number) results. When working with input values where the exponents
are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are
substantially different from the answers provided by quad precision
more than 1% of the time. This error rate may be unacceptable for many
applications that cannot a priori restrict their computations to the
safe range. The proposed method reduces the frequency of
"substantially different" answers by more than 99% for double
precision at a modest cost of performance.
Differences between current gcc methods and the new method will be
described. Then accuracy and performance differences will be discussed.
NOTATION
For all of the following, the notation is:
Input complex values:
a+bi (a= real part, b= imaginary part)
c+di
Output complex value:
e+fi = (a+bi)/(c+di)
For the result tables:
current = current method (SMITH)
b1div = method proposed by Elen Kalda
b2div = alternate method considered by Elen Kalda
new1 = new method using 1 divide and 2 multiplies
new = new method proposed by this patch
DESCRIPTIONS of different complex divide methods:
NAIVE COMPUTATION (-fcx-limited-range):
e = (a*c + b*d)/(c*c + d*d)
f = (b*c - a*d)/(c*c + d*d)
Note that c*c and d*d will overflow or underflow if either
c or d is outside the range 2^-538 to 2^512.
This method is available in gcc when the switch -fcx-limited-range is
used. That switch is also enabled by -ffast-math. Only one who has a
clear understanding of the maximum range of intermediate values
generated by a computation should consider using this switch.
SMITH's METHOD (current libgcc):
if(fabs(c) RBIG) || (FABS (a) > RBIG) || (FABS (b) > RBIG) ) {
a = a * 0.5;
b = b * 0.5;
c = c * 0.5;
d = d * 0.5;
}
/* minimize overflow/underflow issues when c and d are small */
else if (FABS (d) < RMIN2) {
a = a * RMINSCAL;
b = b * RMINSCAL;
c = c * RMINSCAL;
d = d * RMINSCAL;
}
r = c/d; denom = (c*r) + d;
if( r > RMIN ) {
e = (a*r + b) / denom ;
f = (b*r - a) / denom
} else {
e = (c * (a/d) + b) / denom;
f = (c * (b/d) - a) / denom;
}
[ only presenting the fabs(c) < fabs(d) case here, full code in patch. ]
Before any computation of the answer, the co
[PATCH V3] Practical Improvement to libgcc Complex Divide
(Version 3)
(Added in version 3)
Support for half, float, extended, and long double precision has
been added to the prior work for double precision. Since half precision
is computed with float precision as per current libgcc practice,
the enhanced underflow/overflow tests provide no benefit for half
precision and would cost performance. Therefore half precision is
left unchanged.
The existing constants for each precision:
float: FLT_MAX, FLT_MIN;
double: DBL_MAX, DBL_MIN;
extended and/or long double: LDBL_MAX, LDBL_MIN
are used for avoiding the more common overflow/underflow cases.
Additional investigation showed that testing for when both parts of
the denominator had exponents roughly small enough to allow shifting
any subnormal values to normal values, all input values could be
scaled up without risking unnecessary overflow and gaining a clear
improvement in accuracy. The test and scaling values used all fit
within the allowed exponent range for each precision required by the C
standard. The remaining number of troubling results in version 3 is
measurably smaller than in versions 1 and 2.
The timing and precision tables below have been revised appropriately
to match the algorithms used in this version for double precision
and additional tables added to include results for other precisions.
In prior versions, I omitted mention of the bug report that started me
on this project: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=59714
complex division is surprising on targets with FMA (was: on aarch64)
With the proposed method, whether using FMA or not, dividing
1.0+3.0i by 1.0+3.0i correctly returns 1.0+0.0i.
I also have added a reference to Beebe's "The Mathematical Function
Computation Handbook" [4] which was my starting point for research
into better complex divide methods.
(Added for Version 2)
In my initial research, I missed Elen Kalda's proposed patch
https://gcc.gnu.org/legacy-ml/gcc-patches/2019-08/msg01629.html [3]
Thanks to Joseph Myers for providing me with the pointer.
This version includes performance and accuracy comparisions
between Elen's proposed patch and my latest patch version for
double precision.
(from earlier Versions)
The following patch to libgcc/libgcc2.c __divdc3 provides an
opportunity to gain important improvements to the quality of answers
for the default complex divide routine (half, float, double, extended,
long double precisions) when dealing with very large or very small exponents.
The current code correctly implements Smith's method (1962) [1]
further modified by c99's requirements for dealing with NaN (not a
number) results. When working with input values where the exponents
are greater than *_MAX_EXP/2 or less than -(*_MAX_EXP)/2, results are
substantially different from the answers provided by quad precision
more than 1% of the time. This error rate may be unacceptable for many
applications that cannot a priori restrict their computations to the
safe range. The proposed method reduces the frequency of
"substantially different" answers by more than 99% for double
precision at a modest cost of performance.
Differences between current gcc methods and the new method will be
described. Then accuracy and performance differences will be discussed.
NOTATION
For all of the following, the notation is:
Input complex values:
a+bi (a= real part, b= imaginary part)
c+di
Output complex value:
e+fi = (a+bi)/(c+di)
For the result tables:
current = current method (SMITH)
b1div = method proposed by Elen Kalda
b2div = alternate method considered by Elen Kalda
new1 = new method using 1 divide and 2 multiplies
new = new method proposed by this patch
DESCRIPTIONS of different complex divide methods:
NAIVE COMPUTATION (-fcx-limited-range):
e = (a*c + b*d)/(c*c + d*d)
f = (b*c - a*d)/(c*c + d*d)
Note that c*c and d*d will overflow or underflow if either
c or d is outside the range 2^-538 to 2^512.
This method is available in gcc when the switch -fcx-limited-range is
used. That switch is also enabled by -ffast-math. Only one who has a
clear understanding of the maximum range of intermediate values
generated by a computation should consider using this switch.
SMITH's METHOD (current libgcc):
if(fabs(c) RBIG) || (FABS (a) > RBIG) || (FABS (b) > RBIG) ) {
a = a * 0.5;
b = b * 0.5;
c = c * 0.5;
d = d * 0.5;
}
/* minimize overflow/underflow issues when c and d are small */
else if (FABS (d) < RMIN2) {
a = a * RMINSCAL;
b = b * RMINSCAL;
c = c * RMINSCAL;
d = d * RMINSCAL;
}
r = c/d; denom = (c*r) + d;
if( r > RMIN ) {
e = (a*r + b) / denom ;
f = (b*r - a) / denom
} else {
e = (c * (a/d) + b) / denom;
f = (c * (b/d) - a) / denom;
}
[ only presenting the fabs(c) < fabs(d) case here, full code in patch. ]
Before any computation of the answer, the code checks for any input
values near maximum to allow down scaling to avoid overflow. These
scalings almost never harm the accuracy since
Re: [PATCH V2] Practical Improvement to Double Precision Complex Divide
I will study real.c carefully along with the C99 std to determine if I can find useful values for RMIN2 and RMINSCAL for each format which are within range for all instances of that format. A quick skim of real.c shows we have ieee half precision and two arm half precision formats, for example. I appreciate the pointer to real.c as I was sure that information about various platform formats was somewhere in the src tree, but had not found it on my own. - patrick On 6/10/2020 5:39 PM, Joseph Myers wrote: On Wed, 10 Jun 2020, Patrick McGehearty wrote: #ifdef L_divhc3 #define RBIG (correct value for half precision) #define RMIN (correct value for half precision) #define RMIN2 ... (correct value for half precision) #define RMINSCAL ... (correct value for half precision) #endif #ifdef L_divsc3 #define RBIG FLT_MAX #define RMIN FLT_MIN #define RMIN2 ... #define RMINSCAL ... #endif ... would get the job done once I determine (and test) the correct values for all formats. But some of those machine modes can represent several different floating-point formats, depending on the target. There are many more floating-point formats supported in real.c than there are names for floating-point machine modes. Not all the differences are relevant here (those differences relating to default NaNs, for example), but some are.
Re: [PATCH V2] Practical Improvement to Double Precision Complex Divide
Joseph, Thank you again for your prompt and insightful response. It's now clear that I was engaged in wishful thinking that I only needed to handle double precision in my proposed change. I understand now that the text above the routine: - - - - - #if defined(L_divhc3) || defined(L_divsc3) || defined(L_divdc3) \ || defined(L_divxc3) || defined(L_divtc3) - - - - - implies the same code is to be built for L_divhc3 half precision L_divsc3 single precision L_divdc3 double precision L_divxc3 extended precision L_divtc3 quad precision I will need expand my test suites for these other floating point formats and do my best to understand how the libgcc build mechanisms work before I'm ready to resubmit. At first glance, assuming only one of the L_div?c3 names are defined at a time, it seems like something like: #ifdef L_divhc3 #define RBIG (correct value for half precision) #define RMIN (correct value for half precision) #define RMIN2 ... (correct value for half precision) #define RMINSCAL ... (correct value for half precision) #endif #ifdef L_divsc3 #define RBIG FLT_MAX #define RMIN FLT_MIN #define RMIN2 ... #define RMINSCAL ... #endif ... would get the job done once I determine (and test) the correct values for all formats. - Patrick McGehearty On 6/9/2020 6:11 PM, Joseph Myers wrote: On Wed, 10 Jun 2020, Patrick McGehearty wrote: I see your point about other floating point formats. According to the C standard, DOUBLE precision must have a DBL_MAX of at least 1E+37. That is (slightly) greater than 0x1.0p+63. Would #define RMIN2 (0x1.0p-53) #define RMINSCAL (0x1.0p+51) be acceptable? Those will be in range for any architecture that supports the C standard. But this code is used for different machine modes, not just double. In particular, on Arm / AArch64 it may be used to build __divhc3. And those constants are certainly outside the range of HFmode (IEEE binary16). I think you need to work out appropriate logic for what values are correct for a given machine mode, in terms of the parameters of that mode (maximum and minimum exponents etc.). Then, if you can't represent the logic for determining those values (with the correct type, not hardcoding the absence of a suffix which would mean double) directly in C, see how c-cppbuiltin.c predefines extra macros if -fbuilding-libgcc (macros relating to built-in function suffixes and the number of bits in the mantissa of a floating-point mode are included there - it would certainly make sense to add more such macros if necessary).
Re: [PATCH V2] Practical Improvement to Double Precision Complex Divide
A follow up note relating to use of fused multiply add in complex divide:
While reviewing bugs relating to complex divide in libgcc, I
rediscovered bug 59714 - complex division is surprising on targets
with FMA.
The specific concern was when you divide (1.0 + 3.0i) by itself and
fused multiply add was active, you did not get (1.0 + 0.0i) as an
answer. Instead you got 1.0 + epsilon off of zero.
The question raised in that bug report was whether it was appropriate
to use fused multiply add. I.e. on the whole, did it contribute to
greater or less accuracy over a range of cases. No follow up data has
been posted.
My data sets (10M randomly generated values over the full exponent and
mantissa range) showed the following with/without FMA:
Moderate exponent data set (between -512 and 511 for exponent)
error current current new new
exceeds no fma fma no fma fma
1 ulp 0.34492% 0.24715% 0.34492% 0.24715%
2 ulp 0.02360% 0.01764% 0.02360% 0.01764%
Full exponent data set (all exponents generated)
error current current new new
exceeds no fma fma no fma fma
1 ulp 1.90432% 1.87195% 0.23167% 0.19746%
2 ulp 1.70836% 1.70646% 0.04067% 0.03856%
While the differences are not dramatic between using fma and no fma,
they seem fairly clear that fma provides slightly more accurate
results more often than not.
I concur with the current practice of using fma if the architecture
provides it.
As a side note, my new method for complex divide gets the exact result
for dividing (1.0 + 3.0i) by itself whether or not fused multiply add
is used. That's because the calculation of the imaginary part
includes subtracting two nearly equal values which turns out to be
less than DBL_MIN. That means the new algorithm changes the order of
operations which does not involve the use of fused mul-add.
- Patrick McGehearty
On 6/4/2020 6:38 PM, Patrick McGehearty via Gcc-patches wrote:
The following patch to libgcc/libgcc2.c __divdc3 provides an
opportunity to gain important improvements to the quality of answers
for the default double precision complex divide routine when dealing
with very large or very small exponents.
The current code correctly implements Smith's method (1962) [1]
further modified by c99's requirements for dealing with NaN (not a
number) results. When working with input values where the exponents
are greater than 512 (i.e. 2.0^512) or less than -512 (i.e. 2.0^-512),
results are substantially different from the answers provided by quad
precision more than 1% of the time. Since the allowed exponent range
for double precision numbers is -1076 to +1023, the error rate may be
unacceptable for many applications. The proposed method reduces the
frequency of "substantially different" answers by more than 99% at a
modest cost of performance.
Differences between current gcc methods and the new method will be
described. Then accuracy and performance differences will be discussed.
(Added for Version 2)
In my initial research, I missed Elen Kalda's proposed patch
https://gcc.gnu.org/legacy-ml/gcc-patches/2019-08/msg01629.html [3]
Thanks to Joseph Myers for providing me with the pointer.
This version includes performance and accuracy comparisions
between Elen's proposed patch and my latest patch version which
has been modified to eliminate two branches.
NOTATION
For all of the following, the notation is:
Input complex values:
a+bi (a= 64bit real part, b= 64bit imaginary part)
c+di
Output complex value:
e+fi = (a+bi)/(c+di)
For the result tables:
current = current method (SMITH)
b1div = method proposed by Elen Kalda
b2div = alternate method considered by Elen Kalda
new1 = new method using 1 divide and 2 multiplies
new = new method proposed by this patch
DESCRIPTIONS of different complex divide methods:
NAIVE COMPUTATION (-fcx-limited-range):
e = (a*c + b*d)/(c*c + d*d)
f = (b*c - a*d)/(c*c + d*d)
Note that c*c and d*d will overflow or underflow if either
c or d is outside the range 2^-538 to 2^512.
This method is available in gcc when the switch -fcx-limited-range is
used. That switch is also enabled by -ffast-math. Only one who has a
clear understanding of the maximum range of intermediate values
generated by a computation should consider using this switch.
SMITH's METHOD (current libgcc):
if(fabs(c) RBIG) || (FABS (a) > RBIG) || (FABS (b) > RBIG) ) {
a = a * 0.5;
b = b * 0.5;
c = c * 0.5;
d = d * 0.5;
}
/* minimize overflow/underflow issues when c and d are small */
else if (FABS (d) < RMIN2) {
a = a * RMINSCAL;
b = b * RMINSCAL;
c = c * RMINSCAL;
d = d * RMINSCAL;
}
r = c/d; denom = (c*r) + d;
if( r > RMIN ) {
e = (a*r + b) / denom ;
f = (b*r - a) / denom
} else {
e = (c * (a/d) + b) / denom;
f = (c * (b/d) - a) / denom;
}
[ only presenting the fabs(c) < fabs(d) case
Re: [PATCH V2] Practical Improvement to Double Precision Complex Divide
I see your point about other floating point formats.
According to the C standard, DOUBLE precision must
have a DBL_MAX of at least 1E+37. That is (slightly)
greater than 0x1.0p+63.
Would
#define RMIN2 (0x1.0p-53)
#define RMINSCAL (0x1.0p+51)
be acceptable?
Those will be in range for any architecture that supports the C standard.
I estimate using those values will cause a typical complex divide to run
somewhat slower for the full exponent dataset because the test for
whether fabs(d) < RMIN2 [or fabs(c) < RMIN2 ] will be taken more often,
adding four scaling multiplies to the overall operation.
On the positive side, using the new scaling factor significantly reduces
the frequency of 2 ulp or larger inaccuracies for the full range dataset.
The alternative would be to remove the RMIN2 test altogether.
That provides a modest reduction in overhead but increases
the inaccuracy for the 8 ulp case to 0.04147%.
Errors Full Dataset
gtr eq current RMIN2=-512 RMIN2=-53 no RMIN2
==
1 ulp 1.87% 0.19755% 0.16788% 0.22169%
2 ulp 1.70% 0.03857% 0.01007% 0.06182%
8 ulp 1.61% 0.02394% 0.00199% 0.04147%
16 ulp 1.51% 0.01548% 0.00087% 0.02727%
24 ulp 1.41% 0.00945% 0.00042% 0.01694%
52 ulp 1.13% 0.1% 0.1% 0.1%
===
The RMIN2=-53% option greatly reduces 8 ulp or larger errors.
The "No RMIN2" option has about double the error rate
of the initial proposal of RMIN2=-512, but still has a error
rate at the 8 ulp level 40 times lower than the current code.
I welcome any thoughts from the community about the relative importance
of accuracy vs performance tradeoffs at these levels.
- Patrick McGehearty
On 6/4/2020 6:45 PM, Joseph Myers wrote:
On Fri, 5 Jun 2020, Patrick McGehearty wrote:
diff --git a/libgcc/libgcc2.c b/libgcc/libgcc2.c
index e0a9fd7..2a1d3dc 100644
--- a/libgcc/libgcc2.c
+++ b/libgcc/libgcc2.c
@@ -2036,26 +2036,77 @@ CONCAT3(__mul,MODE,3) (MTYPE a, MTYPE b, MTYPE c, MTYPE
d)
CTYPE
CONCAT3(__div,MODE,3) (MTYPE a, MTYPE b, MTYPE c, MTYPE d)
{
+#define RBIG ((DBL_MAX)/2.0)
+#define RMIN (DBL_MIN)
+#define RMIN2 (0x1.0p-512)
+#define RMINSCAL (0x1.0p+510)
This code is used for many different machine modes and floating-point
formats (the type and format corresponding to a particular machine mode
may depend on the target for which GCC is configured). You can't hardcode
particular values specific to DFmode (and to DFmode being IEEE binary64)
here.
[PATCH V2] Practical Improvement to Double Precision Complex Divide
The following patch to libgcc/libgcc2.c __divdc3 provides an
opportunity to gain important improvements to the quality of answers
for the default double precision complex divide routine when dealing
with very large or very small exponents.
The current code correctly implements Smith's method (1962) [1]
further modified by c99's requirements for dealing with NaN (not a
number) results. When working with input values where the exponents
are greater than 512 (i.e. 2.0^512) or less than -512 (i.e. 2.0^-512),
results are substantially different from the answers provided by quad
precision more than 1% of the time. Since the allowed exponent range
for double precision numbers is -1076 to +1023, the error rate may be
unacceptable for many applications. The proposed method reduces the
frequency of "substantially different" answers by more than 99% at a
modest cost of performance.
Differences between current gcc methods and the new method will be
described. Then accuracy and performance differences will be discussed.
(Added for Version 2)
In my initial research, I missed Elen Kalda's proposed patch
https://gcc.gnu.org/legacy-ml/gcc-patches/2019-08/msg01629.html [3]
Thanks to Joseph Myers for providing me with the pointer.
This version includes performance and accuracy comparisions
between Elen's proposed patch and my latest patch version which
has been modified to eliminate two branches.
NOTATION
For all of the following, the notation is:
Input complex values:
a+bi (a= 64bit real part, b= 64bit imaginary part)
c+di
Output complex value:
e+fi = (a+bi)/(c+di)
For the result tables:
current = current method (SMITH)
b1div = method proposed by Elen Kalda
b2div = alternate method considered by Elen Kalda
new1 = new method using 1 divide and 2 multiplies
new = new method proposed by this patch
DESCRIPTIONS of different complex divide methods:
NAIVE COMPUTATION (-fcx-limited-range):
e = (a*c + b*d)/(c*c + d*d)
f = (b*c - a*d)/(c*c + d*d)
Note that c*c and d*d will overflow or underflow if either
c or d is outside the range 2^-538 to 2^512.
This method is available in gcc when the switch -fcx-limited-range is
used. That switch is also enabled by -ffast-math. Only one who has a
clear understanding of the maximum range of intermediate values
generated by a computation should consider using this switch.
SMITH's METHOD (current libgcc):
if(fabs(c) RBIG) || (FABS (a) > RBIG) || (FABS (b) > RBIG) ) {
a = a * 0.5;
b = b * 0.5;
c = c * 0.5;
d = d * 0.5;
}
/* minimize overflow/underflow issues when c and d are small */
else if (FABS (d) < RMIN2) {
a = a * RMINSCAL;
b = b * RMINSCAL;
c = c * RMINSCAL;
d = d * RMINSCAL;
}
r = c/d; denom = (c*r) + d;
if( r > RMIN ) {
e = (a*r + b) / denom ;
f = (b*r - a) / denom
} else {
e = (c * (a/d) + b) / denom;
f = (c * (b/d) - a) / denom;
}
[ only presenting the fabs(c) < fabs(d) case here, full code in patch. ]
Before any computation of the answer, the code checks for near maximum or
near minimum inputs and scale the results to move all values away from
the extremes. If the complex divide can be computed at all without
generating infinities, these scalings will not affect the accuracy
since they are by a power of 2. Values that are over RBIG are
relatively rare but it is easy to test for them and required to avoid
unnecessary overflows.
Testing for RMIN2 reveals when both c and d are less than 2^-512. By
scaling all values by 2^510, the code avoids many underflows in
intermediate computations that otherwise might occur. If scaling a and
b by 2^510 causes either to overflow, then the computation will overflow
whatever method is used.
Next, r (the ratio of c to d) is checked for being near zero. Baudin
and Smith checked r for zero. Checking for values less than DBL_MIN
covers more cases and improves overall accuracy. If r is near zero,
then when it is used in a multiplication, there is a high chance that
the result will underflow to zero, losing significant accuracy. That
underflow can be avoided if the computation is done in a different
order. When r is subnormal, the code replaces a*r (= a*(c/d)) with
((a/d)*c) which is mathematically the same but avoids the unnecessary
underflow.
TEST Data
Two sets of data are presented to test these methods. Both sets
contain 10 million pairs of 64bit complex values. The exponents and
mantissas are generated using multiple calls to random() and then
combining the results. Only values which give results to complex
divide that are representable in 64-bits after being computed in quad
precision are used.
The first data set is labeled "moderate exponents".
The exponent range is limited to -512 to +511.
The second data set is labeled "full exponents".
The exponent range is -1076 to + 1024.
ACCURACY Test results:
Note: All results are based on use of fused multiply-add. If
fused multiply-add is not used, the error rate increases slightly
for
[PATCH] Practical Improvement to Double Precision Complex Divide
The following patch to libgcc/libgcc2.c __divdc3 provides an
opportunity to gain important improvements to the quality of answers
for the default double precision complex divide routine when dealing
with very large or very small exponents.
The current code correctly implements Smith's method (1962) [1]
further modified by c99's requirements for dealing with NaN (not a
number) results. When working with input values where the exponents
are greater than 512 (i.e. 2.0^512) or less than -512 (i.e. 2.0^-512),
results are substantially different from the answers provided by quad
precision more than 1% of the time. Since the allowed exponent range
for double precision numbers is -1076 to +1023, the error rate may be
unacceptable for many applications. The proposed method reduces the
frequency of "substantially different" answers by more than 99% at a
modest cost of performance.
Differences between current gcc methods and the new method will be
described. Then accuracy and performance differences will be discussed.
NOTATION
For all of the following, the notation is:
Input complex values:
a+bi (a= 64bit real part, b= 64bit imaginary part)
c+di
Output complex value:
e+fi = (a+bi)/(c+di)
DESCRIPTIONS of different complex divide methods:
NAIVE COMPUTATION (-fcx-limited-range):
e = (a*c + b*d)/(c*c + d*d)
f = (b*c - a*d)/(c*c + d*d)
Note that c*c and d*d will overflow or underflow if either
c or d is outside the range 2^-538 to 2^512.
This method is available in gcc when the switch -fcx-limited-range is
used. That switch is also enabled by -ffast-math. Only one who has a
clear understanding of the maximum range of intermediate values
generated by a computation should consider using this switch.
SMITH's METHOD (current libgcc):
if(fabs(c) RBIG) || (FABS (a) > RBIG) || (FABS (b) > RBIG) ) {
a = a * 0.5;
b = b * 0.5;
c = c * 0.5;
d = d * 0.5;
}
/* minimize overflow/underflow issues when c and d are small */
else if (FABS (d) < RMIN2) {
a = a * RMINSCAL;
b = b * RMINSCAL;
c = c * RMINSCAL;
d = d * RMINSCAL;
}
r = c/d; denom = (c*r) + d;
if( r > RMIN ) {
e = (a*r + b) / denom ;
f = (b*r - a) / denom
} else {
e = (c * (a/d) + b) / denom;
f = (c * (b/d) - a) / denom;
}
[ only presenting the fabs(c) < fabs(d) case here, full code in patch. ]
Before any computation of the answer, the code checks for near maximum or
near minimum inputs and scale the results to move all values away from
the extremes. If the complex divide can be computed at all without
generating infinities, these scalings will not affect the accuracy
since they are by a power of 2. Values that are over RBIG are
relatively rare but it is easy to test for them and required to avoid
unnecessary overflows.
Testing for RMIN2 reveals when both c and d are less than 2^-512. By
scaling all values by 2^510, the code avoids many underflows in
intermediate computations that otherwise might occur. If scaling a and
b by 2^510 causes either to overflow, then the computation will overflow
whatever method is used.
Next, r (the ratio of c to d) is checked for being near zero. Baudin
and Smith checked r for zero. Checking for values less than DBL_MIN
covers more cases and improves overall accuracy. If r is near zero,
then when it is used in a multiplication, there is a high chance that
the result will underflow to zero, losing significant accuracy. That
underflow can be avoided if the computation is done in a different
order. When r is subnormal, the code replaces a*r (= a*(c/d)) with
((a/d)*c) which is mathematically the same but avoids the unnecessary
underflow.
TEST Data
Two sets of data are presented to test these methods. Both sets
contain 10 million pairs of 64bit complex values. The exponents and
mantissas are generated using multiple calls to random() and then
combining the results. Only values which give results to complex
divide that are representable in 64-bits after being computed in quad
precision are used.
The first data set is labeled "moderate exponents".
The exponent range is limited to -512 to +511.
The second data set is labeled "full exponents".
The exponent range is -1076 to + 1024.
ACCURACY Test results:
Note: All results are based on use of fused multiply-add. If
fused multiply-add is not used, the error rate increases slightly
for the 2 ulp and 8 ulp cases.
The complex divide methods are evaluated by determining what
percentage of values exceed different ulp (units in last place)
levels. If a "2 ulp" test results show 1%, that would mean that 1% of
10,000,000 values (100,000) have either a real or imaginary part that
had a greater than 2 bit difference from the quad precision result.
Results are reported for differences greater than or equal to 2 ulp, 8
ulp, 16 ulp, 24 ulp, and 52 ulp. Even when the patch avoids overflows and
underflows, some input values are expected to have errors due to
normal limitations of floating point
