77#include <complex.h>
88#include <stdio.h>
99
10+ #define RELATIVE_TOLERANCE 1e-9
1011
1112// Returns: the product of a + ib and c + id
1213
@@ -15,6 +16,19 @@ __muldc3(double __a, double __b, double __c, double __d);
1516
1617enum {zero , non_zero , inf , NaN , non_zero_nan };
1718
19+ int check_complex_equal (double _Complex r1 , double _Complex r2 )
20+ {
21+ double max_magnitude = fmax (cabs (r1 ), cabs (r2 ));
22+ double real_diff = fabs (creal (r1 ) - creal (r2 ));
23+ double imag_diff = fabs (cimag (r1 ) - cimag (r2 ));
24+ if (real_diff >= max_magnitude * RELATIVE_TOLERANCE )
25+ return 0 ;
26+ if (imag_diff >= max_magnitude * RELATIVE_TOLERANCE )
27+ return 0 ;
28+
29+ return 1 ;
30+ }
31+
1832int
1933classify (double _Complex x )
2034{
@@ -46,11 +60,15 @@ int test__muldc3(double a, double b, double c, double d)
4660// a, b, c, d, creal(r), cimag(r));
4761 double _Complex dividend ;
4862 double _Complex divisor ;
49-
63+ double _Complex temp ;
64+
5065 __real__ dividend = a ;
5166 __imag__ dividend = b ;
5267 __real__ divisor = c ;
5368 __imag__ divisor = d ;
69+
70+ __real__ temp = a * c - b * d ;
71+ __imag__ temp = a * d + b * c ;
5472
5573 switch (classify (dividend ))
5674 {
@@ -89,7 +107,7 @@ int test__muldc3(double a, double b, double c, double d)
89107 case non_zero :
90108 if (classify (r ) != non_zero )
91109 return 1 ;
92- if (r != a * c - b * d + _Complex_I * ( a * d + b * c ))
110+ if (! check_complex_equal ( r , temp ))
93111 return 1 ;
94112 break ;
95113 case inf :
0 commit comments