lib/testfp.c

   1 /*
   2  *  Test 128-bit floating-point arithmetic on arm64:
   3  *  build with two different compilers and compare the output.
   4  *
   5  *  Copyright (c) 2015 Edmund Grimley Evans
   6  *
   7  * Copying and distribution of this file, with or without modification,
   8  * are permitted in any medium without royalty provided the copyright
   9  * notice and this notice are preserved.  This file is offered as-is,
  10  * without any warranty.
  11  */
  12
  13 #include <stdint.h>
  14 #include <stdio.h>
  15 #include <stdlib.h>
  16 #include <string.h>
  17
  18 #define check(x) ((x) ? (void)0 : check_fail(#x, __FILE__, __LINE__))
  19
  20 void check_fail(const char *assertion, const char *file, unsigned int line)
  21 {
  22     printf("%s:%d: Check (%s) failed.", file, line, assertion);
  23     exit(1);
  24 }
  25
  26 typedef struct {
  27     unsigned long long x0, x1;
  28 } u128_t;
  29
  30 float copy_fi(uint32_t x)
  31 {
  32     float f;
  33     memcpy(&f, &x, 4);
  34     return f;
  35 }
  36
  37 double copy_di(uint64_t x)
  38 {
  39     double f;
  40     memcpy(&f, &x, 8);
  41     return f;
  42 }
  43
  44 long double copy_ldi(u128_t x)
  45 {
  46     long double f;
  47     memcpy(&f, &x, 16);
  48     return f;
  49 }
  50
  51 uint32_t copy_if(float f)
  52 {
  53     uint32_t x;
  54     memcpy(&x, &f, 4);
  55     return x;
  56 }
  57
  58 uint64_t copy_id(double f)
  59 {
  60     uint64_t x;
  61     memcpy(&x, &f, 8);
  62     return x;
  63 }
  64
  65 u128_t copy_ild(long double f)
  66 {
  67     u128_t x;
  68     memcpy(&x, &f, 16);
  69     return x;
  70 }
  71
  72 long double make(int sgn, int exp, uint64_t high, uint64_t low)
  73 {
  74     u128_t x = { low,
  75                  (0x0000ffffffffffff & high) |
  76                  (0x7fff000000000000 & (uint64_t)exp << 48) |
  77                  (0x8000000000000000 & (uint64_t)sgn << 63) };
  78     return copy_ldi(x);
  79 }
  80
  81 void cmp(long double a, long double b)
  82 {
  83     u128_t ax = copy_ild(a);
  84     u128_t bx = copy_ild(b);
  85     int eq = (a == b);
  86     int ne = (a != b);
  87     int lt = (a < b);
  88     int le = (a <= b);
  89     int gt = (a > b);
  90     int ge = (a >= b);
  91
  92     check(eq == 0 || eq == 1);
  93     check(lt == 0 || lt == 1);
  94     check(gt == 0 || gt == 1);
  95     check(ne == !eq && le == (lt | eq) && ge == (gt | eq));
  96     check(eq + lt + gt < 2);
  97
  98     printf("cmp %016llx%016llx %016llx%016llx %d %d %d\n",
  99            ax.x1, ax.x0, bx.x1, bx.x0, lt, eq, gt);
 100 }
 101
 102 void cmps(void)
 103 {
 104     int i, j;
 105
 106     for (i = 0; i < 2; i++)
 107         for (j = 0; j < 2; j++)
 108             cmp(make(i, 0, 0, 0), make(j, 0, 0, 0));
 109
 110     for (i = 0; i < 2; i++) {
 111         for (j = 0; j < 64; j++) {
 112             long double f1 = make(i, 32767, (uint64_t)1 << j, 0);
 113             long double f2 = make(i, 32767, 0, (uint64_t)1 << j);
 114             cmp(f1, 0);
 115             cmp(f2, 0);
 116             cmp(0, f1);
 117             cmp(0, f2);
 118         }
 119     }
 120
 121     for (i = 0; i < 6; i++)
 122         for (j = 0; j < 6; j++)
 123             cmp(make(i & 1, i >> 1, 0, 0),
 124                 make(j & 1, j >> 1, 0, 0));
 125
 126     for (i = 0; i < 2; i++) {
 127         for (j = 0; j < 2; j++) {
 128             int a, b;
 129             for (a = 0; a < 2; a++) {
 130                 for (b = 0; b < 2; b++) {
 131                     cmp(make(i, j, a, b), make(i, j, 0, 0));
 132                     cmp(make(i, j, 0, 0), make(i, j, a, b));
 133                 }
 134             }
 135         }
 136     }
 137 }
 138
 139 void xop(const char *name, long double a, long double b, long double c)
 140 {
 141     u128_t ax = copy_ild(a);
 142     u128_t bx = copy_ild(b);
 143     u128_t cx = copy_ild(c);
 144     printf("%s %016llx%016llx %016llx%016llx %016llx%016llx\n",
 145            name, ax.x1, ax.x0, bx.x1, bx.x0, cx.x1, cx.x0);
 146 }
 147
 148 void fadd(long double a, long double b)
 149 {
 150     xop("add", a, b, a + b);
 151 }
 152
 153 void fsub(long double a, long double b)
 154 {
 155     xop("sub", a, b, a - b);
 156 }
 157
 158 void fmul(long double a, long double b)
 159 {
 160     xop("mul", a, b, a * b);
 161 }
 162
 163 void fdiv(long double a, long double b)
 164 {
 165     xop("div", a, b, a / b);
 166 }
 167
 168 void nanz(void)
 169 {
 170     // Check NaNs:
 171     {
 172         long double x[7];
 173         int i, j, n = 0;
 174         x[n++] = make(0, 32000, 0x95132b76effc, 0xd79035214b4f8d53);
 175         x[n++] = make(1, 32001, 0xbe71d7a51587, 0x30601c6815d6c3ac);
 176         x[n++] = make(0, 32767, 0, 1);
 177         x[n++] = make(0, 32767, (uint64_t)1 << 46, 0);
 178         x[n++] = make(1, 32767, (uint64_t)1 << 47, 0);
 179         x[n++] = make(1, 32767, 0x7596c7099ad5, 0xe25fed2c58f73fc9);
 180         x[n++] = make(0, 32767, 0x835d143360f9, 0x5e315efb35630666);
 181         check(n == sizeof(x) / sizeof(*x));
 182         for (i = 0; i < n; i++) {
 183             for (j = 0; j < n; j++) {
 184                 fadd(x[i], x[j]);
 185                 fsub(x[i], x[j]);
 186                 fmul(x[i], x[j]);
 187                 fdiv(x[i], x[j]);
 188             }
 189         }
 190     }
 191
 192     // Check infinities and zeroes:
 193     {
 194         long double x[6];
 195         int i, j, n = 0;
 196         x[n++] = make(1, 32000, 0x62acda85f700, 0x47b6c9f35edc4044);
 197         x[n++] = make(0, 32001, 0x94b7abf55af7, 0x9f425fe354428e19);
 198         x[n++] = make(0, 32767, 0, 0);
 199         x[n++] = make(1, 32767, 0, 0);
 200         x[n++] = make(0, 0, 0, 0);
 201         x[n++] = make(1, 0, 0, 0);
 202         check(n == sizeof(x) / sizeof(*x));
 203         for (i = 0; i < n; i++) {
 204             for (j = 0; j < n; j++) {
 205                 fadd(x[i], x[j]);
 206                 fsub(x[i], x[j]);
 207                 fmul(x[i], x[j]);
 208                 fdiv(x[i], x[j]);
 209             }
 210         }
 211     }
 212 }
 213
 214 void adds(void)
 215 {
 216     // Check shifting and add/sub:
 217     {
 218         int i;
 219         for (i = -130; i <= 130; i++) {
 220             int s1 = (uint32_t)i % 3 < 1;
 221             int s2 = (uint32_t)i % 5 < 2;
 222             fadd(make(s1, 16384    , 0x502c065e4f71a65d, 0xd2f9bdb031f4f031),
 223                  make(s2, 16384 + i, 0xae267395a9bc1033, 0xb56b5800da1ba448));
 224         }
 225     }
 226
 227     // Check normalisation:
 228     {
 229         uint64_t a0 = 0xc6bab0a6afbef5ed;
 230         uint64_t a1 = 0x4f84136c4a2e9b52;
 231         int ee[] = { 0, 1, 10000 };
 232         int e, i;
 233         for (e = 0; e < sizeof(ee) / sizeof(*ee); e++) {
 234             int exp = ee[e];
 235             fsub(make(0, exp, a1, a0), make(0, 0, 0, 0));
 236             for (i = 63; i >= 0; i--)
 237                 fsub(make(0, exp, a1 | (uint64_t)1 << i >> 1, a0),
 238                      make(0, exp, a1 >> i << i, 0));
 239             for (i = 63; i >=0; i--)
 240                 fsub(make(0, exp, a1, a0 | (uint64_t)1 << i >> 1),
 241                      make(0, exp, a1, a0 >> i << i));
 242         }
 243     }
 244
 245     // Carry/overflow from rounding:
 246     {
 247         fadd(make(0, 114, -1, -1), make(0, 1, 0, 0));
 248         fadd(make(0, 32766, -1, -1), make(0, 32653, 0, 0));
 249         fsub(make(1, 32766, -1, -1), make(0, 32653, 0, 0));
 250     }
 251 }
 252
 253 void muls(void)
 254 {
 255     int i, j;
 256
 257     {
 258         long double max = make(0, 32766, -1, -1);
 259         long double min = make(0, 0, 0, 1);
 260         fmul(max, max);
 261         fmul(max, min);
 262         fmul(min, min);
 263     }
 264
 265     for (i = 117; i > 0; i--)
 266         fmul(make(0, 16268, 0x643dcea76edc, 0xe0877a598403627a),
 267              make(i & 1, i, 0, 0));
 268
 269     fmul(make(0, 16383, -1, -3), make(0, 16383, 0, 1));
 270     // Round to next exponent:
 271     fmul(make(0, 16383, -1, -2), make(0, 16383, 0, 1));
 272     // Round from subnormal to normal:
 273     fmul(make(0, 1, -1, -1), make(0, 16382, 0, 0));
 274
 275     for (i = 0; i < 2; i++)
 276         for (j = 0; j < 112; j++)
 277             fmul(make(0, 16383, (uint64_t)1 << i, 0),
 278                  make(0, 16383,
 279                       j < 64 ? 0 : (uint64_t)1 << (j - 64),
 280                       j < 64 ? (uint64_t)1 << j : 0));
 281 }
 282
 283 void divs(void)
 284 {
 285     int i;
 286
 287     {
 288         long double max = make(0, 32766, -1, -1);
 289         long double min = make(0, 0, 0, 1);
 290         fdiv(max, max);
 291         fdiv(max, min);
 292         fdiv(min, max);
 293         fdiv(min, min);
 294     }
 295
 296     for (i = 0; i < 64; i++)
 297         fdiv(make(0, 16383, -1, -1), make(0, 16383, -1, -(uint64_t)1 << i));
 298     for (i = 0; i < 48; i++)
 299         fdiv(make(0, 16383, -1, -1), make(0, 16383, -(uint64_t)1 << i, 0));
 300 }
 301
 302 void cvtlsw(int32_t a)
 303 {
 304     long double f = a;
 305     u128_t x = copy_ild(f);
 306     printf("cvtlsw %08lx %016llx%016llx\n", (long)(uint32_t)a, x.x1, x.x0);
 307 }
 308
 309 void cvtlsx(int64_t a)
 310 {
 311     long double f = a;
 312     u128_t x = copy_ild(f);
 313     printf("cvtlsx %016llx %016llx%016llx\n",
 314            (long long)(uint64_t)a, x.x1, x.x0);
 315 }
 316
 317 void cvtluw(uint32_t a)
 318 {
 319     long double f = a;
 320     u128_t x = copy_ild(f);
 321     printf("cvtluw %08lx %016llx%016llx\n", (long)a, x.x1, x.x0);
 322 }
 323
 324 void cvtlux(uint64_t a)
 325 {
 326     long double f = a;
 327     u128_t x = copy_ild(f);
 328     printf("cvtlux %016llx %016llx%016llx\n", (long long)a, x.x1, x.x0);
 329 }
 330
 331 void cvtil(long double a)
 332 {
 333     u128_t x = copy_ild(a);
 334     int32_t b1 = a;
 335     int64_t b2 = a;
 336     uint32_t b3 = a;
 337     uint64_t b4 = a;
 338     printf("cvtswl %016llx%016llx %08lx\n",
 339            x.x1, x.x0, (long)(uint32_t)b1);
 340     printf("cvtsxl %016llx%016llx %016llx\n",
 341            x.x1, x.x0, (long long)(uint64_t)b2);
 342     printf("cvtuwl %016llx%016llx %08lx\n",
 343            x.x1, x.x0, (long)b3);
 344     printf("cvtuxl %016llx%016llx %016llx\n",
 345            x.x1, x.x0, (long long)b4);
 346 }
 347
 348 void cvtlf(float a)
 349 {
 350     uint32_t ax = copy_if(a);
 351     long double b = a;
 352     u128_t bx = copy_ild(b);
 353     printf("cvtlf %08lx %016llx%016llx\n", (long)ax, bx.x1, bx.x0);
 354 }
 355
 356 void cvtld(double a)
 357 {
 358     uint64_t ax = copy_id(a);
 359     long double b = a;
 360     u128_t bx = copy_ild(b);
 361     printf("cvtld %016llx %016llx%016llx\n", (long long)ax, bx.x1, bx.x0);
 362 }
 363
 364 void cvtfl(long double a)
 365 {
 366     u128_t ax = copy_ild(a);
 367     float b = a;
 368     uint32_t bx = copy_if(b);
 369     printf("cvtfl %016llx%016llx %08lx\n", ax.x1, ax.x0, (long)bx);
 370 }
 371
 372 void cvtdl(long double a)
 373 {
 374     u128_t ax = copy_ild(a);
 375     double b = a;
 376     uint64_t bx = copy_id(b);
 377     printf("cvtdl %016llx%016llx %016llx\n", ax.x1, ax.x0, (long long)bx);
 378 }
 379
 380 void cvts(void)
 381 {
 382     int i, j;
 383
 384     {
 385         uint32_t x = 0xad040c5b;
 386         cvtlsw(0);
 387         for (i = 0; i < 31; i++)
 388             cvtlsw(x >> (31 - i));
 389         for (i = 0; i < 31; i++)
 390             cvtlsw(-(x >> (31 - i)));
 391         cvtlsw(0x80000000);
 392     }
 393     {
 394         uint64_t x = 0xb630a248cad9afd2;
 395         cvtlsx(0);
 396         for (i = 0; i < 63; i++)
 397             cvtlsx(x >> (63 - i));
 398         for (i = 0; i < 63; i++)
 399             cvtlsx(-(x >> (63 - i)));
 400         cvtlsx(0x8000000000000000);
 401     }
 402     {
 403         uint32_t x = 0xad040c5b;
 404         cvtluw(0);
 405         for (i = 0; i < 32; i++)
 406             cvtluw(x >> (31 - i));
 407     }
 408     {
 409         uint64_t x = 0xb630a248cad9afd2;
 410         cvtlux(0);
 411         for (i = 0; i < 64; i++)
 412             cvtlux(x >> (63 - i));
 413     }
 414
 415     for (i = 0; i < 2; i++) {
 416         cvtil(make(i, 32767, 0, 1));
 417         cvtil(make(i, 32767, (uint64_t)1 << 47, 0));
 418         cvtil(make(i, 32767, 123, 456));
 419         cvtil(make(i, 32767, 0, 0));
 420         cvtil(make(i, 16382, -1, -1));
 421         cvtil(make(i, 16383, -1, -1));
 422         cvtil(make(i, 16384, 0x7fffffffffff, -1));
 423         cvtil(make(i, 16384, 0x800000000000, 0));
 424         for (j = 0; j < 68; j++)
 425             cvtil(make(i, 16381 + j, 0xd4822c0a10ec, 0x1fe2f8b2669f5c9d));
 426     }
 427
 428     cvtlf(copy_fi(0x00000000));
 429     cvtlf(copy_fi(0x456789ab));
 430     cvtlf(copy_fi(0x7f800000));
 431     cvtlf(copy_fi(0x7f923456));
 432     cvtlf(copy_fi(0x7fdbcdef));
 433     cvtlf(copy_fi(0x80000000));
 434     cvtlf(copy_fi(0xabcdef12));
 435     cvtlf(copy_fi(0xff800000));
 436     cvtlf(copy_fi(0xff923456));
 437     cvtlf(copy_fi(0xffdbcdef));
 438
 439     cvtld(copy_di(0x0000000000000000));
 440     cvtld(copy_di(0x456789abcdef0123));
 441     cvtld(copy_di(0x7ff0000000000000));
 442     cvtld(copy_di(0x7ff123456789abcd));
 443     cvtld(copy_di(0x7ffabcdef1234567));
 444     cvtld(copy_di(0x8000000000000000));
 445     cvtld(copy_di(0xcdef123456789abc));
 446     cvtld(copy_di(0xfff0000000000000));
 447     cvtld(copy_di(0xfff123456789abcd));
 448     cvtld(copy_di(0xfffabcdef1234567));
 449
 450     for (i = 0; i < 2; i++) {                   \
 451         cvtfl(make(i, 0, 0, 0));
 452         cvtfl(make(i, 16232, -1, -1));
 453         cvtfl(make(i, 16233, 0, 0));
 454         cvtfl(make(i, 16233, 0, 1));
 455         cvtfl(make(i, 16383, 0xab0ffd000000, 0));
 456         cvtfl(make(i, 16383, 0xab0ffd000001, 0));
 457         cvtfl(make(i, 16383, 0xab0ffeffffff, 0));
 458         cvtfl(make(i, 16383, 0xab0fff000000, 0));
 459         cvtfl(make(i, 16383, 0xab0fff000001, 0));
 460         cvtfl(make(i, 16510, 0xfffffeffffff, -1));
 461         cvtfl(make(i, 16510, 0xffffff000000, 0));
 462         cvtfl(make(i, 16511, 0, 0));
 463         cvtfl(make(i, 32767, 0, 0));
 464         cvtfl(make(i, 32767, 0, 1));
 465         cvtfl(make(i, 32767, 0x4cbe01ac5f40, 0x75cee3c6afbb00b5));
 466         cvtfl(make(i, 32767, 0x800000000000, 1));
 467         cvtfl(make(i, 32767, 0xa11caaaf6a52, 0x696033e871eab099));
 468     }
 469
 470     for (i = 0; i < 2; i++) {
 471         cvtdl(make(i, 0, 0, 0));
 472         cvtdl(make(i, 15307, -1, -1));
 473         cvtdl(make(i, 15308, 0, 0));
 474         cvtdl(make(i, 15308, 0, 1));
 475         cvtdl(make(i, 16383, 0xabc123abc0ff, 0xe800000000000000));
 476         cvtdl(make(i, 16383, 0xabc123abc0ff, 0xe800000000000001));
 477         cvtdl(make(i, 16383, 0xabc123abc0ff, 0xf7ffffffffffffff));
 478         cvtdl(make(i, 16383, 0xabc123abc0ff, 0xf800000000000000));
 479         cvtdl(make(i, 16383, 0xabc123abc0ff, 0xf800000000000001));
 480         cvtdl(make(i, 17406, 0xffffffffffff, 0xf7ffffffffffffff));
 481         cvtdl(make(i, 17406, 0xffffffffffff, 0xf800000000000000));
 482         cvtdl(make(i, 17407, 0, 0));
 483         cvtdl(make(i, 32767, 0, 0));
 484         cvtdl(make(i, 32767, 0, 1));
 485         cvtdl(make(i, 32767, 0x4cbe01ac5f40, 0x75cee3c6afbb00b5));
 486         cvtdl(make(i, 32767, 0x800000000000, 1));
 487         cvtdl(make(i, 32767, 0xa11caaaf6a52, 0x696033e871eab099));
 488     }
 489 }
 490
 491 void tests(void)
 492 {
 493     cmps();
 494     nanz();
 495     adds();
 496     muls();
 497     divs();
 498     cvts();
 499 }
 500
 501 int main()
 502 {
 503 #ifdef __aarch64__
 504     tests();
 505 #else
 506     printf("This test program is intended for a little-endian architecture\n"
 507            "with an IEEE-standard 128-bit long double.\n");
 508 #endif
 509     return 0;
 510 }