1 // Copyright Naoki Shibata and contributors 2010 - 2020.
2 // Distributed under the Boost Software License, Version 1.0.
3 // (See accompanying file LICENSE.txt or copy at
4 // http://www.boost.org/LICENSE_1_0.txt)
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <stdint.h>
9 #include <mpfr.h>
10 #include <time.h>
11 #include <float.h>
12 #include <limits.h>
13 #include <assert.h>
14
15 #include <math.h>
16 #include <quadmath.h>
17
18 #define _GNU_SOURCE
19 #include <unistd.h>
20 #include <sys/syscall.h>
21 #include <linux/random.h>
22
23 #include "sleef.h"
24
25 #include "f128util.h"
26
27 #define DORENAME
28 #include "rename.h"
29
30 #define POSITIVE_INFINITY INFINITY
31 #define NEGATIVE_INFINITY (-INFINITY)
32
isnumberq(Sleef_quad x)33 int isnumberq(Sleef_quad x) { return !isinfq(x) && !isnanq(x); }
isPlusZeroq(Sleef_quad x)34 int isPlusZeroq(Sleef_quad x) { return x == 0 && copysignq(1, x) == 1; }
isMinusZeroq(Sleef_quad x)35 int isMinusZeroq(Sleef_quad x) { return x == 0 && copysignq(1, x) == -1; }
36
37 mpfr_t fra, frb, frc, frd;
38
countULP(Sleef_quad d,mpfr_t c)39 double countULP(Sleef_quad d, mpfr_t c) {
40 Sleef_quad c2 = mpfr_get_f128(c, GMP_RNDN);
41 if (c2 == 0 && d != 0) return 10000;
42 //if (isPlusZeroq(c2) && !isPlusZeroq(d)) return 10003;
43 //if (isMinusZeroq(c2) && !isMinusZeroq(d)) return 10004;
44 if (isnanq(c2) && isnanq(d)) return 0;
45 if (isnanq(c2) || isnanq(d)) return 10001;
46 if (c2 == POSITIVE_INFINITY && d == POSITIVE_INFINITY) return 0;
47 if (c2 == NEGATIVE_INFINITY && d == NEGATIVE_INFINITY) return 0;
48 if (!isnumberq(c2) && !isnumberq(d)) return 0;
49
50 int e;
51 frexpq(mpfr_get_f128(c, GMP_RNDN), &e);
52 mpfr_set_f128(frb, fmaxq(ldexpq(1.0, e-113), FLT128_DENORM_MIN), GMP_RNDN);
53
54 mpfr_set_f128(frd, d, GMP_RNDN);
55 mpfr_sub(fra, frd, c, GMP_RNDN);
56 mpfr_div(fra, fra, frb, GMP_RNDN);
57 double u = fabs(mpfr_get_d(fra, GMP_RNDN));
58
59 return u;
60 }
61
countULP2(Sleef_quad d,mpfr_t c)62 double countULP2(Sleef_quad d, mpfr_t c) {
63 Sleef_quad c2 = mpfr_get_f128(c, GMP_RNDN);
64 if (c2 == 0 && d != 0) return 10000;
65 //if (isPlusZeroq(c2) && !isPlusZeroq(d)) return 10003;
66 //if (isMinusZeroq(c2) && !isMinusZeroq(d)) return 10004;
67 if (isnanq(c2) && isnanq(d)) return 0;
68 if (isnanq(c2) || isnanq(d)) return 10001;
69 if (c2 == POSITIVE_INFINITY && d == POSITIVE_INFINITY) return 0;
70 if (c2 == NEGATIVE_INFINITY && d == NEGATIVE_INFINITY) return 0;
71 if (!isnumberq(c2) && !isnumberq(d)) return 0;
72
73 int e;
74 frexpq(mpfr_get_f128(c, GMP_RNDN), &e);
75 mpfr_set_f128(frb, fmaxq(ldexpq(1.0, e-113), FLT128_MIN), GMP_RNDN);
76
77 mpfr_set_f128(frd, d, GMP_RNDN);
78 mpfr_sub(fra, frd, c, GMP_RNDN);
79 mpfr_div(fra, fra, frb, GMP_RNDN);
80 double u = fabs(mpfr_get_d(fra, GMP_RNDN));
81
82 return u;
83 }
84
85 typedef union {
86 Sleef_quad d;
87 __int128 u128;
88 uint64_t u[2];
89 } conv_t;
90
rnd()91 Sleef_quad rnd() {
92 conv_t c;
93 switch(random() & 15) {
94 case 0: return INFINITY;
95 case 1: return -INFINITY;
96 }
97 syscall(SYS_getrandom, &c.u128, sizeof(c.u128), 0);
98 return c.d;
99 }
100
rnd_fr()101 Sleef_quad rnd_fr() {
102 conv_t c;
103 do {
104 syscall(SYS_getrandom, &c.u128, sizeof(c.u128), 0);
105 } while(!isnumberq(c.d));
106 return c.d;
107 }
108
rnd_zo()109 Sleef_quad rnd_zo() {
110 conv_t c;
111 do {
112 syscall(SYS_getrandom, &c.u128, sizeof(c.u128), 0);
113 } while(!isnumberq(c.d) || c.d < -1 || 1 < c.d);
114 return c.d;
115 }
116
sinpifr(mpfr_t ret,Sleef_quad d)117 void sinpifr(mpfr_t ret, Sleef_quad d) {
118 mpfr_t frpi, frd;
119 mpfr_inits(frpi, frd, NULL);
120
121 mpfr_const_pi(frpi, GMP_RNDN);
122 mpfr_set_d(frd, 1.0, GMP_RNDN);
123 mpfr_mul(frpi, frpi, frd, GMP_RNDN);
124 mpfr_set_f128(frd, d, GMP_RNDN);
125 mpfr_mul(frd, frpi, frd, GMP_RNDN);
126 mpfr_sin(ret, frd, GMP_RNDN);
127
128 mpfr_clears(frpi, frd, NULL);
129 }
130
cospifr(mpfr_t ret,Sleef_quad d)131 void cospifr(mpfr_t ret, Sleef_quad d) {
132 mpfr_t frpi, frd;
133 mpfr_inits(frpi, frd, NULL);
134
135 mpfr_const_pi(frpi, GMP_RNDN);
136 mpfr_set_d(frd, 1.0, GMP_RNDN);
137 mpfr_mul(frpi, frpi, frd, GMP_RNDN);
138 mpfr_set_f128(frd, d, GMP_RNDN);
139 mpfr_mul(frd, frpi, frd, GMP_RNDN);
140 mpfr_cos(ret, frd, GMP_RNDN);
141
142 mpfr_clears(frpi, frd, NULL);
143 }
144
main(int argc,char ** argv)145 int main(int argc,char **argv)
146 {
147 mpfr_t frw, frx, fry, frz;
148
149 mpfr_set_default_prec(2048);
150 mpfr_inits(fra, frb, frc, frd, frw, frx, fry, frz, NULL);
151
152 conv_t cd;
153 Sleef_quad d, t, d2, zo;
154
155 int cnt, ecnt = 0;
156
157 srandom(time(NULL));
158
159 #if 0
160 cd.d = M_PIq;
161 mpfr_set_f128(frx, cd.d, GMP_RNDN);
162 cd.u128 += 3;
163 printf("%g\n", countULP2(cd.d, frx));
164 #endif
165
166 const Sleef_quad rangemax = 1e+9;
167
168 for(cnt = 0;ecnt < 1000;cnt++) {
169 switch(cnt & 7) {
170 case 0:
171 d = rnd();
172 d2 = rnd();
173 zo = rnd();
174 break;
175 case 1:
176 cd.d = rint((2 * (double)random() / RAND_MAX - 1) * 1e+10) * M_PI_4;
177 cd.u128 += (random() & 0xff) - 0x7f;
178 d = cd.d;
179 d2 = rnd();
180 zo = rnd();
181 break;
182 default:
183 d = rnd_fr();
184 d2 = rnd_fr();
185 zo = rnd_zo();
186 break;
187 }
188
189 Sleef_quad2 sc = xsincospiq_u05(d);
190 Sleef_quad2 sc2 = xsincospiq_u35(d);
191
192 {
193 const double rangemax2 = 1e+9;
194
195 sinpifr(frx, d);
196
197 double u0 = countULP2(t = sc.x, frx);
198
199 if (u0 != 0 && ((fabs(d) <= rangemax2 && u0 > 0.505) || fabs(t) > 1 || !isnumberq(t))) {
200 printf("Pure C sincospiq_u05 sin arg="); printf128(d); printf(" ulp=%.20g\n", u0);
201 fflush(stdout); ecnt++;
202 }
203
204 double u1 = countULP2(t = sc2.x, frx);
205
206 if (u1 != 0 && ((fabs(d) <= rangemax2 && u1 > 2.0) || fabs(t) > 1 || !isnumberq(t))) {
207 printf("Pure C sincospiq_u35 sin arg=%.30Lg ulp=%.20g\n", (long double)d, u1);
208 fflush(stdout); ecnt++;
209 }
210
211 }
212
213 {
214 const double rangemax2 = 1e+9;
215
216 cospifr(frx, d);
217
218 double u0 = countULP2(t = sc.y, frx);
219
220 if (u0 != 0 && ((fabs(d) <= rangemax2 && u0 > 0.505) || fabs(t) > 1 || !isnumberq(t))) {
221 printf("Pure C sincospiq_u05 cos arg=%.30Lg ulp=%.20g\n", (long double)d, u0);
222 fflush(stdout); ecnt++;
223 }
224
225 double u1 = countULP2(t = sc.y, frx);
226
227 if (u1 != 0 && ((fabs(d) <= rangemax2 && u1 > 2.0) || fabs(t) > 1 || !isnumberq(t))) {
228 printf("Pure C sincospiq_u35 cos arg=%.30Lg ulp=%.20g\n", (long double)d, u1);
229 fflush(stdout); ecnt++;
230 }
231
232 }
233
234 #if 0
235 double2 sc = xsincos(d);
236 double2 sc2 = xsincos_u1(d);
237
238 {
239 mpfr_set_d(frx, d, GMP_RNDN);
240 mpfr_sin(frx, frx, GMP_RNDN);
241
242 double u0 = countULP(t = xsin(d), frx);
243
244 if ((fabs(d) <= rangemax && u0 > 3.5) || fabs(t) > 1 || !isnumberq(t)) {
245 printf("Pure C sin arg=%.20g ulp=%.20g\n", d, u0);
246 fflush(stdout); ecnt++;
247 }
248
249 double u1 = countULP(sc.x, frx);
250
251 if ((fabs(d) <= rangemax && u1 > 3.5) || fabs(t) > 1 || !isnumberq(t)) {
252 printf("Pure C sincos sin arg=%.20g ulp=%.20g\n", d, u1);
253 fflush(stdout); ecnt++;
254 }
255
256 double u2 = countULP(t = xsin_u1(d), frx);
257
258 if ((fabs(d) <= rangemax && u2 > 1) || fabs(t) > 1 || !isnumberq(t)) {
259 printf("Pure C sin_u1 arg=%.20g ulp=%.20g\n", d, u2);
260 fflush(stdout); ecnt++;
261 }
262
263 double u3 = countULP(t = sc2.x, frx);
264
265 if ((fabs(d) <= rangemax && u3 > 1) || fabs(t) > 1 || !isnumberq(t)) {
266 printf("Pure C sincos_u1 sin arg=%.20g ulp=%.20g\n", d, u3);
267 fflush(stdout); ecnt++;
268 }
269 }
270
271 {
272 mpfr_set_d(frx, d, GMP_RNDN);
273 mpfr_cos(frx, frx, GMP_RNDN);
274
275 double u0 = countULP(t = xcos(d), frx);
276
277 if ((fabs(d) <= rangemax && u0 > 3.5) || fabs(t) > 1 || !isnumberq(t)) {
278 printf("Pure C cos arg=%.20g ulp=%.20g\n", d, u0);
279 fflush(stdout); ecnt++;
280 }
281
282 double u1 = countULP(t = sc.y, frx);
283
284 if ((fabs(d) <= rangemax && u1 > 3.5) || fabs(t) > 1 || !isnumberq(t)) {
285 printf("Pure C sincos cos arg=%.20g ulp=%.20g\n", d, u1);
286 fflush(stdout); ecnt++;
287 }
288
289 double u2 = countULP(t = xcos_u1(d), frx);
290
291 if ((fabs(d) <= rangemax && u2 > 1) || fabs(t) > 1 || !isnumberq(t)) {
292 printf("Pure C cos_u1 arg=%.20g ulp=%.20g\n", d, u2);
293 fflush(stdout); ecnt++;
294 }
295
296 double u3 = countULP(t = sc2.y, frx);
297
298 if ((fabs(d) <= rangemax && u3 > 1) || fabs(t) > 1 || !isnumberq(t)) {
299 printf("Pure C sincos_u1 cos arg=%.20g ulp=%.20g\n", d, u3);
300 fflush(stdout); ecnt++;
301 }
302 }
303
304 {
305 mpfr_set_d(frx, d, GMP_RNDN);
306 mpfr_tan(frx, frx, GMP_RNDN);
307
308 double u0 = countULP(t = xtan(d), frx);
309
310 if ((fabs(d) < 1e+7 && u0 > 3.5) || (fabs(d) <= rangemax && u0 > 5) || isnan(t)) {
311 printf("Pure C tan arg=%.20g ulp=%.20g\n", d, u0);
312 fflush(stdout); ecnt++;
313 }
314
315 double u1 = countULP(t = xtan_u1(d), frx);
316
317 if ((fabs(d) <= rangemax && u1 > 1) || isnan(t)) {
318 printf("Pure C tan_u1 arg=%.20g ulp=%.20g\n", d, u1);
319 fflush(stdout); ecnt++;
320 }
321 }
322
323 d = rnd_fr();
324 double d2 = rnd_fr(), zo = rnd_zo();
325
326 {
327 mpfr_set_d(frx, fabs(d), GMP_RNDN);
328 mpfr_log(frx, frx, GMP_RNDN);
329
330 double u0 = countULP(t = xlog(fabs(d)), frx);
331
332 if (u0 > 3.5) {
333 printf("Pure C log arg=%.20g ulp=%.20g\n", d, u0);
334 fflush(stdout); ecnt++;
335 }
336
337 double u1 = countULP(t = xlog_u1(fabs(d)), frx);
338
339 if (u1 > 1) {
340 printf("Pure C log_u1 arg=%.20g ulp=%.20g\n", d, u1);
341 fflush(stdout); ecnt++;
342 }
343 }
344
345 {
346 mpfr_set_d(frx, fabs(d), GMP_RNDN);
347 mpfr_log10(frx, frx, GMP_RNDN);
348
349 double u0 = countULP(t = xlog10(fabs(d)), frx);
350
351 if (u0 > 1) {
352 printf("Pure C log10 arg=%.20g ulp=%.20g\n", d, u0);
353 fflush(stdout); ecnt++;
354 }
355 }
356
357 {
358 mpfr_set_d(frx, d, GMP_RNDN);
359 mpfr_log1p(frx, frx, GMP_RNDN);
360
361 double u0 = countULP(t = xlog1p(d), frx);
362
363 if ((-1 <= d && d <= 1e+307 && u0 > 1) ||
364 (d < -1 && !isnan(t)) ||
365 (d > 1e+307 && !(u0 <= 1 || isinf(t)))) {
366 printf("Pure C log1p arg=%.20g ulp=%.20g\n", d, u0);
367 printf("%g\n", t);
368 fflush(stdout); ecnt++;
369 }
370 }
371
372 {
373 mpfr_set_d(frx, d, GMP_RNDN);
374 mpfr_exp(frx, frx, GMP_RNDN);
375
376 double u0 = countULP(t = xexp(d), frx);
377
378 if (u0 > 1) {
379 printf("Pure C exp arg=%.20g ulp=%.20g\n", d, u0);
380 fflush(stdout); ecnt++;
381 }
382 }
383
384 {
385 mpfr_set_d(frx, d, GMP_RNDN);
386 mpfr_exp2(frx, frx, GMP_RNDN);
387
388 double u0 = countULP(t = xexp2(d), frx);
389
390 if (u0 > 1) {
391 printf("Pure C exp2 arg=%.20g ulp=%.20g\n", d, u0);
392 fflush(stdout); ecnt++;
393 }
394 }
395
396 {
397 mpfr_set_d(frx, d, GMP_RNDN);
398 mpfr_exp10(frx, frx, GMP_RNDN);
399
400 double u0 = countULP(t = xexp10(d), frx);
401
402 if (u0 > 1) {
403 printf("Pure C exp10 arg=%.20g ulp=%.20g\n", d, u0);
404 fflush(stdout); ecnt++;
405 }
406 }
407
408 {
409 mpfr_set_d(frx, d, GMP_RNDN);
410 mpfr_expm1(frx, frx, GMP_RNDN);
411
412 double u0 = countULP(t = xexpm1(d), frx);
413
414 if (u0 > 1) {
415 printf("Pure C expm1 arg=%.20g ulp=%.20g\n", d, u0);
416 fflush(stdout); ecnt++;
417 }
418 }
419
420 {
421 mpfr_set_d(frx, d, GMP_RNDN);
422 mpfr_set_d(fry, d2, GMP_RNDN);
423 mpfr_pow(frx, fry, frx, GMP_RNDN);
424
425 double u0 = countULP(t = xpow(d2, d), frx);
426
427 if (u0 > 1) {
428 printf("Pure C pow arg=%.20g, %.20g ulp=%.20g\n", d2, d, u0);
429 fflush(stdout); ecnt++;
430 }
431 }
432
433 {
434 mpfr_set_d(frx, d, GMP_RNDN);
435 mpfr_cbrt(frx, frx, GMP_RNDN);
436
437 double u0 = countULP(t = xcbrt(d), frx);
438
439 if (u0 > 3.5) {
440 printf("Pure C cbrt arg=%.20g ulp=%.20g\n", d, u0);
441 fflush(stdout); ecnt++;
442 }
443
444 double u1 = countULP(t = xcbrt_u1(d), frx);
445
446 if (u1 > 1) {
447 printf("Pure C cbrt_u1 arg=%.20g ulp=%.20g\n", d, u1);
448 fflush(stdout); ecnt++;
449 }
450 }
451
452 {
453 mpfr_set_d(frx, zo, GMP_RNDN);
454 mpfr_asin(frx, frx, GMP_RNDN);
455
456 double u0 = countULP(t = xasin(zo), frx);
457
458 if (u0 > 3.5) {
459 printf("Pure C asin arg=%.20g ulp=%.20g\n", d, u0);
460 fflush(stdout); ecnt++;
461 }
462
463 double u1 = countULP(t = xasin_u1(zo), frx);
464
465 if (u1 > 1) {
466 printf("Pure C asin_u1 arg=%.20g ulp=%.20g\n", d, u1);
467 fflush(stdout); ecnt++;
468 }
469 }
470
471 {
472 mpfr_set_d(frx, zo, GMP_RNDN);
473 mpfr_acos(frx, frx, GMP_RNDN);
474
475 double u0 = countULP(t = xacos(zo), frx);
476
477 if (u0 > 3.5) {
478 printf("Pure C acos arg=%.20g ulp=%.20g\n", d, u0);
479 fflush(stdout); ecnt++;
480 }
481
482 double u1 = countULP(t = xacos_u1(zo), frx);
483
484 if (u1 > 1) {
485 printf("Pure C acos_u1 arg=%.20g ulp=%.20g\n", d, u1);
486 fflush(stdout); ecnt++;
487 }
488 }
489
490 {
491 mpfr_set_d(frx, d, GMP_RNDN);
492 mpfr_atan(frx, frx, GMP_RNDN);
493
494 double u0 = countULP(t = xatan(d), frx);
495
496 if (u0 > 3.5) {
497 printf("Pure C atan arg=%.20g ulp=%.20g\n", d, u0);
498 fflush(stdout); ecnt++;
499 }
500
501 double u1 = countULP(t = xatan_u1(d), frx);
502
503 if (u1 > 1) {
504 printf("Pure C atan_u1 arg=%.20g ulp=%.20g\n", d, u1);
505 fflush(stdout); ecnt++;
506 }
507 }
508
509 {
510 mpfr_set_d(frx, d, GMP_RNDN);
511 mpfr_set_d(fry, d2, GMP_RNDN);
512 mpfr_atan2(frx, fry, frx, GMP_RNDN);
513
514 double u0 = countULP(t = xatan2(d2, d), frx);
515
516 if (u0 > 3.5) {
517 printf("Pure C atan2 arg=%.20g, %.20g ulp=%.20g\n", d2, d, u0);
518 fflush(stdout); ecnt++;
519 }
520
521 double u1 = countULP2(t = xatan2_u1(d2, d), frx);
522
523 if (u1 > 1) {
524 printf("Pure C atan2_u1 arg=%.20g, %.20g ulp=%.20g\n", d2, d, u1);
525 fflush(stdout); ecnt++;
526 }
527 }
528
529 {
530 mpfr_set_d(frx, d, GMP_RNDN);
531 mpfr_sinh(frx, frx, GMP_RNDN);
532
533 double u0 = countULP(t = xsinh(d), frx);
534
535 if ((fabs(d) <= 709 && u0 > 1) ||
536 (d > 709 && !(u0 <= 1 || (isinf(t) && t > 0))) ||
537 (d < -709 && !(u0 <= 1 || (isinf(t) && t < 0)))) {
538 printf("Pure C sinh arg=%.20g ulp=%.20g\n", d, u0);
539 fflush(stdout); ecnt++;
540 }
541 }
542
543 {
544 mpfr_set_d(frx, d, GMP_RNDN);
545 mpfr_cosh(frx, frx, GMP_RNDN);
546
547 double u0 = countULP(t = xcosh(d), frx);
548
549 if ((fabs(d) <= 709 && u0 > 1) || !(u0 <= 1 || (isinf(t) && t > 0))) {
550 printf("Pure C cosh arg=%.20g ulp=%.20g\n", d, u0);
551 fflush(stdout); ecnt++;
552 }
553 }
554
555 {
556 mpfr_set_d(frx, d, GMP_RNDN);
557 mpfr_tanh(frx, frx, GMP_RNDN);
558
559 double u0 = countULP(t = xtanh(d), frx);
560
561 if (u0 > 1) {
562 printf("Pure C tanh arg=%.20g ulp=%.20g\n", d, u0);
563 fflush(stdout); ecnt++;
564 }
565 }
566
567 {
568 mpfr_set_d(frx, d, GMP_RNDN);
569 mpfr_asinh(frx, frx, GMP_RNDN);
570
571 double u0 = countULP(t = xasinh(d), frx);
572
573 if ((fabs(d) < sqrt(DBL_MAX) && u0 > 1) ||
574 (d >= sqrt(DBL_MAX) && !(u0 <= 1 || (isinf(t) && t > 0))) ||
575 (d <= -sqrt(DBL_MAX) && !(u0 <= 1 || (isinf(t) && t < 0)))) {
576 printf("Pure C asinh arg=%.20g ulp=%.20g\n", d, u0);
577 fflush(stdout); ecnt++;
578 }
579 }
580
581 {
582 mpfr_set_d(frx, d, GMP_RNDN);
583 mpfr_acosh(frx, frx, GMP_RNDN);
584
585 double u0 = countULP(t = xacosh(d), frx);
586
587 if ((fabs(d) < sqrt(DBL_MAX) && u0 > 1) ||
588 (d >= sqrt(DBL_MAX) && !(u0 <= 1 || (isinf(t) && t > 0))) ||
589 (d <= -sqrt(DBL_MAX) && !isnan(t))) {
590 printf("Pure C acosh arg=%.20g ulp=%.20g\n", d, u0);
591 printf("%.20g\n", t);
592 fflush(stdout); ecnt++;
593 }
594 }
595
596 {
597 mpfr_set_d(frx, d, GMP_RNDN);
598 mpfr_atanh(frx, frx, GMP_RNDN);
599
600 double u0 = countULP(t = xatanh(d), frx);
601
602 if (u0 > 1) {
603 printf("Pure C atanh arg=%.20g ulp=%.20g\n", d, u0);
604 fflush(stdout); ecnt++;
605 }
606 }
607 #endif
608 }
609 }
610