[builtins] Fix overflow issue for complex division with big numbers

[zsrkmyn]

Summary

__div*c3 functions get overflow with large numbers, e.g. (DBL_MAX + DBL_MAXi) / (0.5DBL_MAX + 0.5DBL_MAXi) for __divdc3.
This patch tries to fix them.

[zsrkmyn]

ping

[craig.topper]

@scanon, is this something you can help with?

[scanon]

lib/builtins/divdc3.c

34

We don't have a clear definition of the required semantics for these builtins, but should FMA formation be allowed here?

As written, it looks like we're at the mercy of the fp-contract flags passed while building compiler-rt, which means that the exact cancellation assumed by the tests won't pass if the default is ever "on". I would suggest either relaxing the tests or adding an explicit #pragma STDC FP_CONTRACT OFF in these functions.

The latter is the simpler option, and doesn't have any real downside; we can discuss relaxing it later if we want.

test/builtins/Unit/divdc3_test.c

368

It would be good to build up these tests a bit more. In particular for double, we should include the test cases from Baudin & Smith's "A Robust Complex Division in Scilab" (https://arxiv.org/abs/1210.4539) for double; there are only 10, but they cover most of the hardest-to-get-right finite-result cases for complex division.

They are trying to implement a division with a *componentwise* error-bound. That's a non-goal for compiler-rt, so we shouldn't expect to get them exactly right, but we should be able to get them with small relative error.

[scanon]

lib/builtins/divdc3.c

33

This rescaling strategy is still prone to overflow. Consider a = b = DBL_MAX, c = d = 1. The rescaled a*c + b*d will still overflow, because max(c,d) is 1, but the result should be finite.

In general, I don't think that you can implement a rescaled division like this without considering the scaling of both (a,b) and (c,d); it doesn't suffice to only use the scaling of (c,d) alone.

lib/builtins/divsc3.c

21

On targets where we have hardware double-precision, it would be a huge win for float to simply promote to double and skip rescaling entirely. We still need the inf/nan checks, but for finite results this always gives the right answer.

This can be a follow-on patch, however.

lib/builtins/divdc3.c

 // Returns: the quotient of (a + ib) / (c + id)
 
 COMPILER_RT_ABI Dcomplex __divdc3(double __a, double __b, double __c,
                                   double __d) {
   int __ilogbw = 0;
   double __logbw = __compiler_rt_logb(crt_fmax(crt_fabs(__c), crt_fabs(__d)));
   if (crt_isfinite(__logbw)) {
     __ilogbw = (int)__logbw;
+    __a = crt_scalbn(__a, -__ilogbw);
+    __b = crt_scalbn(__b, -__ilogbw);
     __c = crt_scalbn(__c, -__ilogbw);
     __d = crt_scalbn(__d, -__ilogbw);
   }
   double __denom = __c * __c + __d * __d;
   Dcomplex z;
-  COMPLEX_REAL(z) = crt_scalbn((__a * __c + __b * __d) / __denom, -__ilogbw);
-  COMPLEX_IMAGINARY(z) =
-      crt_scalbn((__b * __c - __a * __d) / __denom, -__ilogbw);
+  COMPLEX_REAL(z) = (__a * __c + __b * __d) / __denom;
+  COMPLEX_IMAGINARY(z) = (__b * __c - __a * __d) / __denom;
   if (crt_isnan(COMPLEX_REAL(z)) && crt_isnan(COMPLEX_IMAGINARY(z))) {
     if ((__denom == 0.0) && (!crt_isnan(__a) || !crt_isnan(__b))) {
       COMPLEX_REAL(z) = crt_copysign(CRT_INFINITY, __c) * __a;
       COMPLEX_IMAGINARY(z) = crt_copysign(CRT_INFINITY, __c) * __b;
     } else if ((crt_isinf(__a) || crt_isinf(__b)) && crt_isfinite(__c) &&
                crt_isfinite(__d)) {
       __a = crt_copysign(crt_isinf(__a) ? 1.0 : 0.0, __a);
       __b = crt_copysign(crt_isinf(__b) ? 1.0 : 0.0, __b);

lib/builtins/divsc3.c

 #define SINGLE_PRECISION
 #include "fp_lib.h"
 #include "int_lib.h"
 #include "int_math.h"
 
 // Returns: the quotient of (a + ib) / (c + id)
 
 COMPILER_RT_ABI Fcomplex __divsc3(float __a, float __b, float __c, float __d) {
   int __ilogbw = 0;
   float __logbw =
       __compiler_rt_logbf(crt_fmaxf(crt_fabsf(__c), crt_fabsf(__d)));
   if (crt_isfinite(__logbw)) {
     __ilogbw = (int)__logbw;
+    __a = crt_scalbnf(__a, -__ilogbw);
+    __b = crt_scalbnf(__b, -__ilogbw);
     __c = crt_scalbnf(__c, -__ilogbw);
     __d = crt_scalbnf(__d, -__ilogbw);
   }
   float __denom = __c * __c + __d * __d;
   Fcomplex z;
-  COMPLEX_REAL(z) = crt_scalbnf((__a * __c + __b * __d) / __denom, -__ilogbw);
-  COMPLEX_IMAGINARY(z) =
-      crt_scalbnf((__b * __c - __a * __d) / __denom, -__ilogbw);
+  COMPLEX_REAL(z) = (__a * __c + __b * __d) / __denom;
+  COMPLEX_IMAGINARY(z) = (__b * __c - __a * __d) / __denom;
   if (crt_isnan(COMPLEX_REAL(z)) && crt_isnan(COMPLEX_IMAGINARY(z))) {
     if ((__denom == 0) && (!crt_isnan(__a) || !crt_isnan(__b))) {
       COMPLEX_REAL(z) = crt_copysignf(CRT_INFINITY, __c) * __a;
       COMPLEX_IMAGINARY(z) = crt_copysignf(CRT_INFINITY, __c) * __b;
     } else if ((crt_isinf(__a) || crt_isinf(__b)) && crt_isfinite(__c) &&
                crt_isfinite(__d)) {
       __a = crt_copysignf(crt_isinf(__a) ? 1 : 0, __a);
       __b = crt_copysignf(crt_isinf(__b) ? 1 : 0, __b);

lib/builtins/divtc3.c

 
 COMPILER_RT_ABI Lcomplex __divtc3(long double __a, long double __b,
                                   long double __c, long double __d) {
   int __ilogbw = 0;
   long double __logbw =
       __compiler_rt_logbl(crt_fmaxl(crt_fabsl(__c), crt_fabsl(__d)));
   if (crt_isfinite(__logbw)) {
     __ilogbw = (int)__logbw;
+    __a = crt_scalbnl(__a, -__ilogbw);
+    __b = crt_scalbnl(__b, -__ilogbw);
     __c = crt_scalbnl(__c, -__ilogbw);
     __d = crt_scalbnl(__d, -__ilogbw);
   }
   long double __denom = __c * __c + __d * __d;
   Lcomplex z;
-  COMPLEX_REAL(z) = crt_scalbnl((__a * __c + __b * __d) / __denom, -__ilogbw);
-  COMPLEX_IMAGINARY(z) =
-      crt_scalbnl((__b * __c - __a * __d) / __denom, -__ilogbw);
+  COMPLEX_REAL(z) = (__a * __c + __b * __d) / __denom;
+  COMPLEX_IMAGINARY(z) = (__b * __c - __a * __d) / __denom;
   if (crt_isnan(COMPLEX_REAL(z)) && crt_isnan(COMPLEX_IMAGINARY(z))) {
     if ((__denom == 0.0) && (!crt_isnan(__a) || !crt_isnan(__b))) {
       COMPLEX_REAL(z) = crt_copysignl(CRT_INFINITY, __c) * __a;
       COMPLEX_IMAGINARY(z) = crt_copysignl(CRT_INFINITY, __c) * __b;
     } else if ((crt_isinf(__a) || crt_isinf(__b)) && crt_isfinite(__c) &&
                crt_isfinite(__d)) {
       __a = crt_copysignl(crt_isinf(__a) ? 1.0 : 0.0, __a);
       __b = crt_copysignl(crt_isinf(__b) ? 1.0 : 0.0, __b);

lib/builtins/divxc3.c

 // Returns: the quotient of (a + ib) / (c + id)
 
 COMPILER_RT_ABI Lcomplex __divxc3(long double __a, long double __b,
                                   long double __c, long double __d) {
   int __ilogbw = 0;
   long double __logbw = crt_logbl(crt_fmaxl(crt_fabsl(__c), crt_fabsl(__d)));
   if (crt_isfinite(__logbw)) {
     __ilogbw = (int)__logbw;
+    __a = crt_scalbnl(__a, -__ilogbw);
+    __b = crt_scalbnl(__b, -__ilogbw);
     __c = crt_scalbnl(__c, -__ilogbw);
     __d = crt_scalbnl(__d, -__ilogbw);
   }
   long double __denom = __c * __c + __d * __d;
   Lcomplex z;
-  COMPLEX_REAL(z) = crt_scalbnl((__a * __c + __b * __d) / __denom, -__ilogbw);
-  COMPLEX_IMAGINARY(z) =
-      crt_scalbnl((__b * __c - __a * __d) / __denom, -__ilogbw);
+  COMPLEX_REAL(z) = (__a * __c + __b * __d) / __denom;
+  COMPLEX_IMAGINARY(z) = (__b * __c - __a * __d) / __denom;
   if (crt_isnan(COMPLEX_REAL(z)) && crt_isnan(COMPLEX_IMAGINARY(z))) {
     if ((__denom == 0) && (!crt_isnan(__a) || !crt_isnan(__b))) {
       COMPLEX_REAL(z) = crt_copysignl(CRT_INFINITY, __c) * __a;
       COMPLEX_IMAGINARY(z) = crt_copysignl(CRT_INFINITY, __c) * __b;
     } else if ((crt_isinf(__a) || crt_isinf(__b)) && crt_isfinite(__c) &&
                crt_isfinite(__d)) {
       __a = crt_copysignl(crt_isinf(__a) ? 1 : 0, __a);
       __b = crt_copysignl(crt_isinf(__b) ? 1 : 0, __b);

test/builtins/Unit/divdc3_test.c

classify(double _Complex x)
     return non_zero;
 }
 
 int test__divdc3(double a, double b, double c, double d)
 {
     double _Complex r = __divdc3(a, b, c, d);
 //     printf("test__divdc3(%f, %f, %f, %f) = %f + I%f\n",
 //             a, b, c, d, creal(r), cimag(r));
-  double _Complex dividend;
-  double _Complex divisor;
-  
-  __real__ dividend = a;
-  __imag__ dividend = b;
-  __real__ divisor = c;
-  __imag__ divisor = d;
-  
+    double _Complex dividend;
+    double _Complex divisor;
+    __real__ dividend = a;
+    __imag__ dividend = b;
+    __real__ divisor = c;
+    __imag__ divisor = d;
     switch (classify(dividend))
     {
     case zero:
         switch (classify(divisor))
         {
         case zero:
             if (classify(r) != NaN)
                 return 1;
case non_zero_nan:
             break;
         case non_zero_nan:
             if (classify(r) != NaN)
                 return 1;
             break;
         }
         break;
     }
-    
     return 0;
 }
 
 double x[][2] =
 {
     { 1.e-6,  1.e-6},
     {-1.e-6,  1.e-6},
     {-1.e-6, -1.e-6},
double x[][2] =
     {+0., INFINITY},
     {0.5, INFINITY},
     {1, INFINITY},
     {2, INFINITY},
     {INFINITY, INFINITY}
 
 };
 
+int test_bignum__divdc3(double a, double b, double c, double d, double _Complex r)
+{
+    return __divdc3(a, b, c, d) != r;
+}
+
+double x_bignum[][5] =
+{
+    {__DBL_MAX__, __DBL_MAX__, __DBL_MAX__/2, __DBL_MAX__/2, 2.},
+    {-__DBL_MAX__, -__DBL_MAX__, __DBL_MAX__/2, __DBL_MAX__/2, -2.},
+};
+
 int main()
 {
-    const unsigned N = sizeof(x) / sizeof(x[0]);
+    unsigned N = sizeof(x) / sizeof(x[0]);
     unsigned i, j;
     for (i = 0; i < N; ++i)
     {
         for (j = 0; j < N; ++j)
         {
             if (test__divdc3(x[i][0], x[i][1], x[j][0], x[j][1]))
                 return 1;
         }
     }
 
+    N = sizeof(x_bignum) / sizeof(x_bignum[0]);
+    for (i = 0; i < N; ++i)
+    {
+        double *t = x_bignum[i];
+        if (test_bignum__divdc3(t[0], t[1], t[2], t[3], t[4]))
+            return 1;
+    }
+
     return 0;
 }

test/builtins/Unit/divsc3_test.c

classify(float _Complex x)
     return non_zero;
 }
 
 int test__divsc3(float a, float b, float c, float d)
 {
     float _Complex r = __divsc3(a, b, c, d);
 //     printf("test__divsc3(%f, %f, %f, %f) = %f + I%f\n",
 //             a, b, c, d, crealf(r), cimagf(r));
-  float _Complex dividend;
-  float _Complex divisor;
-  
-  __real__ dividend = a;
-  __imag__ dividend = b;
-  __real__ divisor = c;
-  __imag__ divisor = d;
-  
+    float _Complex dividend;
+    float _Complex divisor;
+    __real__ dividend = a;
+    __imag__ dividend = b;
+    __real__ divisor = c;
+    __imag__ divisor = d;
     switch (classify(dividend))
     {
     case zero:
         switch (classify(divisor))
         {
         case zero:
             if (classify(r) != NaN)
                 return 1;
case non_zero_nan:
             break;
         case non_zero_nan:
             if (classify(r) != NaN)
                 return 1;
             break;
         }
         break;
     }
-    
     return 0;
 }
 
 float x[][2] =
 {
     { 1.e-6,  1.e-6},
     {-1.e-6,  1.e-6},
     {-1.e-6, -1.e-6},
float x[][2] =
     {+0., INFINITY},
     {0.5, INFINITY},
     {1, INFINITY},
     {2, INFINITY},
     {INFINITY, INFINITY}
 
 };
 
+int test_bignum__divsc3(float a, float b, float c, float d, float _Complex r)
+{
+    return __divsc3(a, b, c, d) != r;
+}
+
+float x_bignum[][5] =
+{
+    {__FLT_MAX__, __FLT_MAX__, __FLT_MAX__/2, __FLT_MAX__/2, 2.},
+    {-__FLT_MAX__, -__FLT_MAX__, __FLT_MAX__/2, __FLT_MAX__/2, -2.},
+};
+
 int main()
 {
-    const unsigned N = sizeof(x) / sizeof(x[0]);
+    unsigned N = sizeof(x) / sizeof(x[0]);
     unsigned i, j;
     for (i = 0; i < N; ++i)
     {
         for (j = 0; j < N; ++j)
         {
             if (test__divsc3(x[i][0], x[i][1], x[j][0], x[j][1]))
                 return 1;
         }
     }
 
+    N = sizeof(x_bignum) / sizeof(x_bignum[0]);
+    for (i = 0; i < N; ++i)
+    {
+        float *t = x_bignum[i];
+        if (test_bignum__divsc3(t[0], t[1], t[2], t[3], t[4]))
+            return 1;
+    }
+
     return 0;
 }

test/builtins/Unit/divtc3_test.c

classify(long double _Complex x)
     return non_zero;
 }
 
 int test__divtc3(long double a, long double b, long double c, long double d)
 {
     long double _Complex r = __divtc3(a, b, c, d);
 //      printf("test__divtc3(%Lf, %Lf, %Lf, %Lf) = %Lf + I%Lf\n",
 //              a, b, c, d, creall(r), cimagl(r));
-  
-  long double _Complex dividend;
-  long double _Complex divisor;
-  
-  __real__ dividend = a;
-  __imag__ dividend = b;
-  __real__ divisor = c;
-  __imag__ divisor = d;
-  
+    long double _Complex dividend;
+    long double _Complex divisor;
+    __real__ dividend = a;
+    __imag__ dividend = b;
+    __real__ divisor = c;
+    __imag__ divisor = d;
     switch (classify(dividend))
     {
     case zero:
         switch (classify(divisor))
         {
         case zero:
             if (classify(r) != NaN)
                 return 1;
case non_zero_nan:
             break;
         case non_zero_nan:
             if (classify(r) != NaN)
                 return 1;
             break;
         }
         break;
     }
-    
     return 0;
 }
 
 long double x[][2] =
 {
     { 1.e-6,  1.e-6},
     {-1.e-6,  1.e-6},
     {-1.e-6, -1.e-6},
long double x[][2] =
     {+0., INFINITY},
     {0.5, INFINITY},
     {1, INFINITY},
     {2, INFINITY},
     {INFINITY, INFINITY}
 
 };
 
+int test_bignum__divtc3(long double a, long double b, long double c, long double d, long double _Complex r)
+{
+    return __divtc3(a, b, c, d) != r;
+}
+
+long double x_bignum[][5] =
+{
+    {__LDBL_MAX__, __LDBL_MAX__, __LDBL_MAX__/2, __LDBL_MAX__/2, 2.},
+    {-__LDBL_MAX__, -__LDBL_MAX__, __LDBL_MAX__/2, __LDBL_MAX__/2, -2.},
+};
+
 int main()
 {
-    const unsigned N = sizeof(x) / sizeof(x[0]);
+    unsigned N = sizeof(x) / sizeof(x[0]);
     unsigned i, j;
     for (i = 0; i < N; ++i)
     {
         for (j = 0; j < N; ++j)
         {
             if (test__divtc3(x[i][0], x[i][1], x[j][0], x[j][1]))
                 return 1;
         }
     }
 
-//  printf("No errors found.\n");
+    N = sizeof(x_bignum) / sizeof(x_bignum[0]);
+    for (i = 0; i < N; ++i)
+    {
+        long double *t = x_bignum[i];
+        if (test_bignum__divtc3(t[0], t[1], t[2], t[3], t[4]))
+            return 1;
+    }
+
+//    printf("No errors found.\n");
     return 0;
 }

test/builtins/Unit/divxc3_test.c

classify(long double _Complex x)
     return non_zero;
 }
 
 int test__divxc3(long double a, long double b, long double c, long double d)
 {
     long double _Complex r = __divxc3(a, b, c, d);
 //     printf("test__divxc3(%Lf, %Lf, %Lf, %Lf) = %Lf + I%Lf\n",
 //             a, b, c, d, creall(r), cimagl(r));
-  long double _Complex dividend;
-  long double _Complex divisor;
-  
-  __real__ dividend = a;
-  __imag__ dividend = b;
-  __real__ divisor = c;
-  __imag__ divisor = d;
-  
+    long double _Complex dividend;
+    long double _Complex divisor;
+    __real__ dividend = a;
+    __imag__ dividend = b;
+    __real__ divisor = c;
+    __imag__ divisor = d;
     switch (classify(dividend))
     {
     case zero:
         switch (classify(divisor))
         {
         case zero:
             if (classify(r) != NaN)
                 return 1;
case non_zero_nan:
             break;
         case non_zero_nan:
             if (classify(r) != NaN)
                 return 1;
             break;
         }
         break;
     }
-    
     return 0;
 }
 
 long double x[][2] =
 {
     { 1.e-6,  1.e-6},
     {-1.e-6,  1.e-6},
     {-1.e-6, -1.e-6},
long double x[][2] =
     {+0., INFINITY},
     {0.5, INFINITY},
     {1, INFINITY},
     {2, INFINITY},
     {INFINITY, INFINITY}
 
 };
 
+int test_bignum__divxc3(long double a, long double b, long double c, long double d, long double _Complex r)
+{
+    return __divxc3(a, b, c, d) != r;
+}
+
+long double x_bignum[][5] =
+{
+    {__LDBL_MAX__, __LDBL_MAX__, __LDBL_MAX__/2, __LDBL_MAX__/2, 2.},
+    {-__LDBL_MAX__, -__LDBL_MAX__, __LDBL_MAX__/2, __LDBL_MAX__/2, -2.},
+};
+
 #endif
 
 int main()
 {
 #if !_ARCH_PPC
-    const unsigned N = sizeof(x) / sizeof(x[0]);
+    unsigned N = sizeof(x) / sizeof(x[0]);
     unsigned i, j;
     for (i = 0; i < N; ++i)
     {
         for (j = 0; j < N; ++j)
         {
             if (test__divxc3(x[i][0], x[i][1], x[j][0], x[j][1]))
                 return 1;
         }
     }
 
+    N = sizeof(x_bignum) / sizeof(x_bignum[0]);
+    for (i = 0; i < N; ++i)
+    {
+        long double *t = x_bignum[i];
+        if (test_bignum__divxc3(t[0], t[1], t[2], t[3], t[4]))
+            return 1;
+    }
+
 #else
     printf("skipped\n");
 #endif
     return 0;
 }