7
7
#include <complex.h>
8
8
#include <stdio.h>
9
9
10
+ #define RELATIVE_TOLERANCE 1e-9
10
11
11
12
// Returns: the product of a + ib and c + id
12
13
@@ -15,6 +16,19 @@ __muldc3(double __a, double __b, double __c, double __d);
15
16
16
17
enum {zero , non_zero , inf , NaN , non_zero_nan };
17
18
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
+
18
32
int
19
33
classify (double _Complex x )
20
34
{
@@ -46,11 +60,15 @@ int test__muldc3(double a, double b, double c, double d)
46
60
// a, b, c, d, creal(r), cimag(r));
47
61
double _Complex dividend ;
48
62
double _Complex divisor ;
49
-
63
+ double _Complex temp ;
64
+
50
65
__real__ dividend = a ;
51
66
__imag__ dividend = b ;
52
67
__real__ divisor = c ;
53
68
__imag__ divisor = d ;
69
+
70
+ __real__ temp = a * c - b * d ;
71
+ __imag__ temp = a * d + b * c ;
54
72
55
73
switch (classify (dividend ))
56
74
{
@@ -89,7 +107,7 @@ int test__muldc3(double a, double b, double c, double d)
89
107
case non_zero :
90
108
if (classify (r ) != non_zero )
91
109
return 1 ;
92
- if (r != a * c - b * d + _Complex_I * ( a * d + b * c ))
110
+ if (! check_complex_equal ( r , temp ))
93
111
return 1 ;
94
112
break ;
95
113
case inf :
0 commit comments