4
4
#include "DD.h"
5
5
6
6
long double __gcc_qadd (long double x , long double y );
7
+ int memcmp (const void * , const void * , __typeof__ (sizeof (0 )));
7
8
8
9
double testAccuracy ();
9
10
int testEdgeCases ();
10
11
11
12
int main (int argc , char * argv []) {
12
13
if (testEdgeCases ())
13
14
return 1 ;
14
-
15
+
15
16
if (testAccuracy () > 1.0 )
16
17
return 1 ;
17
-
18
+
18
19
return 0 ;
19
20
}
20
21
@@ -114,11 +115,11 @@ int testEdgeCases() {
114
115
b .lo = edgeCases [i ].ylo ;
115
116
r .hi = edgeCases [i ].rhi ;
116
117
r .lo = edgeCases [i ].rlo ;
117
-
118
+
118
119
int error ;
119
-
120
+
120
121
DD c = { .ld = __gcc_qadd (a .ld , b .ld ) };
121
-
122
+
122
123
if (r .hi != r .hi ) {
123
124
if (c .hi == c .hi )
124
125
error = 1 ;
@@ -130,134 +131,134 @@ int testEdgeCases() {
130
131
131
132
else if (r .hi == 0.0 )
132
133
error = memcmp (& c , & r , sizeof (DD ));
133
-
134
+
134
135
else
135
136
error = ((c .hi != r .hi ) || (c .lo != r .lo ));
136
-
137
+
137
138
if (error ) {
138
139
printf ("Error on edge case %a + %a: expected (%a, %a), got (%a, %a).\n" , a .hi , b .hi , r .hi , r .lo , c .hi , c .lo );
139
140
return 1 ;
140
141
}
141
142
}
142
-
143
+
143
144
return 0 ;
144
145
}
145
146
146
147
147
148
/*
148
149
149
150
Code for generating the test cases, requires the mpfr package to run.
150
-
151
+
151
152
#include <stdio.h>
152
153
#include <stdlib.h>
153
154
#include <mpfr.h>
154
-
155
+
155
156
#ifdef __x86_64__
156
157
#define randlength 2
157
158
#else
158
159
#define randlength 4
159
160
#endif
160
-
161
-
161
+
162
+
162
163
int main(int argc, char *argv[]) {
163
-
164
+
164
165
MPFR_DECL_INIT(a, 106);
165
166
MPFR_DECL_INIT(b, 106);
166
167
MPFR_DECL_INIT(c, 106);
167
-
168
+
168
169
MPFR_DECL_INIT(tmp, 53);
169
-
170
+
170
171
int exponent_range = atoi(argv[1]);
171
-
172
+
172
173
int i;
173
174
for (i=0; i<128; ++i) {
174
175
mpfr_random2(a, randlength, exponent_range);
175
176
mpfr_random2(b, randlength, exponent_range);
176
177
mpfr_add(c, a, b, GMP_RNDN);
177
-
178
+
178
179
double ahi = mpfr_get_d(a, GMP_RNDN);
179
180
mpfr_set_d(tmp, ahi, GMP_RNDN);
180
181
mpfr_sub(tmp, a, tmp, GMP_RNDN);
181
182
double alo = mpfr_get_d(tmp, GMP_RNDN);
182
183
printf("{%0.13a, %0.13a, ", ahi, alo);
183
-
184
+
184
185
double bhi = mpfr_get_d(b, GMP_RNDN);
185
186
mpfr_set_d(tmp, bhi, GMP_RNDN);
186
187
mpfr_sub(tmp, b, tmp, GMP_RNDN);
187
188
double blo = mpfr_get_d(tmp, GMP_RNDN);
188
189
printf("%0.13a, %0.13a, ", bhi, blo);
189
-
190
+
190
191
double chi = mpfr_get_d(c, GMP_RNDN);
191
192
mpfr_set_d(tmp, chi, GMP_RNDN);
192
193
mpfr_sub(tmp, c, tmp, GMP_RNDN);
193
194
double clo = mpfr_get_d(tmp, GMP_RNDN);
194
195
printf("%0.13a, %0.13a},\n", chi, clo);
195
-
196
+
196
197
mpfr_neg(b, b, GMP_RNDN);
197
198
mpfr_add(c, a, b, GMP_RNDN);
198
-
199
+
199
200
ahi = mpfr_get_d(a, GMP_RNDN);
200
201
mpfr_set_d(tmp, ahi, GMP_RNDN);
201
202
mpfr_sub(tmp, a, tmp, GMP_RNDN);
202
203
alo = mpfr_get_d(tmp, GMP_RNDN);
203
204
printf("{%0.13a, %0.13a, ", ahi, alo);
204
-
205
+
205
206
bhi = mpfr_get_d(b, GMP_RNDN);
206
207
mpfr_set_d(tmp, bhi, GMP_RNDN);
207
208
mpfr_sub(tmp, b, tmp, GMP_RNDN);
208
209
blo = mpfr_get_d(tmp, GMP_RNDN);
209
210
printf("%0.13a, %0.13a, ", bhi, blo);
210
-
211
+
211
212
chi = mpfr_get_d(c, GMP_RNDN);
212
213
mpfr_set_d(tmp, chi, GMP_RNDN);
213
214
mpfr_sub(tmp, c, tmp, GMP_RNDN);
214
215
clo = mpfr_get_d(tmp, GMP_RNDN);
215
216
printf("%0.13a, %0.13a},\n", chi, clo);
216
-
217
+
217
218
mpfr_neg(a, a, GMP_RNDN);
218
219
mpfr_neg(b, b, GMP_RNDN);
219
220
mpfr_add(c, a, b, GMP_RNDN);
220
-
221
+
221
222
ahi = mpfr_get_d(a, GMP_RNDN);
222
223
mpfr_set_d(tmp, ahi, GMP_RNDN);
223
224
mpfr_sub(tmp, a, tmp, GMP_RNDN);
224
225
alo = mpfr_get_d(tmp, GMP_RNDN);
225
226
printf("{%0.13a, %0.13a, ", ahi, alo);
226
-
227
+
227
228
bhi = mpfr_get_d(b, GMP_RNDN);
228
229
mpfr_set_d(tmp, bhi, GMP_RNDN);
229
230
mpfr_sub(tmp, b, tmp, GMP_RNDN);
230
231
blo = mpfr_get_d(tmp, GMP_RNDN);
231
232
printf("%0.13a, %0.13a, ", bhi, blo);
232
-
233
+
233
234
chi = mpfr_get_d(c, GMP_RNDN);
234
235
mpfr_set_d(tmp, chi, GMP_RNDN);
235
236
mpfr_sub(tmp, c, tmp, GMP_RNDN);
236
237
clo = mpfr_get_d(tmp, GMP_RNDN);
237
238
printf("%0.13a, %0.13a},\n", chi, clo);
238
-
239
+
239
240
mpfr_neg(b, b, GMP_RNDN);
240
241
mpfr_add(c, a, b, GMP_RNDN);
241
-
242
+
242
243
ahi = mpfr_get_d(a, GMP_RNDN);
243
244
mpfr_set_d(tmp, ahi, GMP_RNDN);
244
245
mpfr_sub(tmp, a, tmp, GMP_RNDN);
245
246
alo = mpfr_get_d(tmp, GMP_RNDN);
246
247
printf("{%0.13a, %0.13a, ", ahi, alo);
247
-
248
+
248
249
bhi = mpfr_get_d(b, GMP_RNDN);
249
250
mpfr_set_d(tmp, bhi, GMP_RNDN);
250
251
mpfr_sub(tmp, b, tmp, GMP_RNDN);
251
252
blo = mpfr_get_d(tmp, GMP_RNDN);
252
253
printf("%0.13a, %0.13a, ", bhi, blo);
253
-
254
+
254
255
chi = mpfr_get_d(c, GMP_RNDN);
255
256
mpfr_set_d(tmp, chi, GMP_RNDN);
256
257
mpfr_sub(tmp, c, tmp, GMP_RNDN);
257
258
clo = mpfr_get_d(tmp, GMP_RNDN);
258
259
printf("%0.13a, %0.13a},\n", chi, clo);
259
260
}
260
-
261
+
261
262
return 0;
262
263
}
263
264
@@ -1813,27 +1814,27 @@ const int numAccuracyTests = sizeof(accuracyTests) / sizeof(struct testVector);
1813
1814
double testAccuracy () {
1814
1815
int i ;
1815
1816
DD a , b , c , r ;
1816
-
1817
+
1817
1818
double worstUlps = 0.5 ;
1818
-
1819
+
1819
1820
for (i = 0 ; i < numAccuracyTests ; ++ i ) {
1820
1821
a .hi = accuracyTests [i ].xhi ;
1821
1822
a .lo = accuracyTests [i ].xlo ;
1822
1823
b .hi = accuracyTests [i ].yhi ;
1823
1824
b .lo = accuracyTests [i ].ylo ;
1824
1825
r .hi = accuracyTests [i ].rhi ;
1825
1826
r .lo = accuracyTests [i ].rlo ;
1826
-
1827
+
1827
1828
DD c = { .ld = __gcc_qadd (a .ld , b .ld ) };
1828
-
1829
+
1829
1830
double error = __builtin_fabs (((r .hi - c .hi ) + r .lo ) - c .lo );
1830
-
1831
+
1831
1832
if (error != 0.0 ) {
1832
-
1833
+
1833
1834
int exponent = ilogb (r .hi );
1834
1835
exponent = (exponent < -1022 ? -1022 : exponent );
1835
1836
double ulpError = scalbn (error , 106 - exponent );
1836
-
1837
+
1837
1838
if (ulpError > worstUlps ) {
1838
1839
#ifdef PRINT_ACCURACY_INFORMATION
1839
1840
printf ("New worst rounding error for (%a,%a) + (%a,%a):\n" , a .hi , a .lo , b .hi , b .lo );
@@ -1845,7 +1846,7 @@ double testAccuracy() {
1845
1846
}
1846
1847
}
1847
1848
}
1848
-
1849
+
1849
1850
return worstUlps ;
1850
1851
}
1851
1852
0 commit comments