contrib/ntp/tests/libntp/lfpfunc.c

   1 #include "config.h"
   2
   3 #include "ntp_stdlib.h"
   4 #include "ntp_fp.h"
   5
   6 #include "unity.h"
   7
   8 #include <float.h>
   9 #include <math.h>
  10
  11
  12 /*
  13    replaced:    TEST_ASSERT_EQUAL_MEMORY(&a, &b, sizeof(a))
  14    with:        TEST_ASSERT_EQUAL_l_fp(a, b).
  15    It's safer this way, because structs can be compared even if they
  16    aren't initiated with memset (due to padding bytes).
  17 */
  18 #define TEST_ASSERT_EQUAL_l_fp(a, b) {                                  \
  19         TEST_ASSERT_EQUAL_MESSAGE(a.l_i, b.l_i, "Field l_i");           \
  20         TEST_ASSERT_EQUAL_UINT_MESSAGE(a.l_uf, b.l_uf, "Field l_uf");   \
  21 }
  22
  23
  24 typedef int bool; // typedef enum { FALSE, TRUE } boolean; -> can't use this because TRUE and FALSE are already defined
  25
  26
  27 typedef struct  {
  28         uint32_t h, l;
  29 } lfp_hl;
  30
  31
  32 int     l_fp_scmp(const l_fp first, const l_fp second);
  33 int     l_fp_ucmp(const l_fp first, l_fp second);
  34 l_fp    l_fp_init(int32 i, u_int32 f);
  35 l_fp    l_fp_add(const l_fp first, const l_fp second);
  36 l_fp    l_fp_subtract(const l_fp first, const l_fp second);
  37 l_fp    l_fp_negate(const l_fp first);
  38 l_fp    l_fp_abs(const l_fp first);
  39 int     l_fp_signum(const l_fp first);
  40 double  l_fp_convert_to_double(const l_fp first);
  41 l_fp    l_fp_init_from_double( double rhs);
  42 void    l_fp_swap(l_fp * first, l_fp *second);
  43 bool    l_isgt(const l_fp first, const l_fp second);
  44 bool    l_isgtu(const l_fp first, const l_fp second);
  45 bool    l_ishis(const l_fp first, const l_fp second);
  46 bool    l_isgeq(const l_fp first, const l_fp second);
  47 bool    l_isequ(const l_fp first, const l_fp second);
  48 double  eps(double d);
  49
  50
  51 void test_AdditionLR(void);
  52 void test_AdditionRL(void);
  53 void test_SubtractionLR(void);
  54 void test_SubtractionRL(void);
  55 void test_Negation(void);
  56 void test_Absolute(void);
  57 void test_FDF_RoundTrip(void);
  58 void test_SignedRelOps(void);
  59 void test_UnsignedRelOps(void);
  60
  61
  62 static int cmp_work(u_int32 a[3], u_int32 b[3]);
  63
  64 //----------------------------------------------------------------------
  65 // reference comparision
  66 // This is implementad as a full signed MP-subtract in 3 limbs, where
  67 // the operands are zero or sign extended before the subtraction is
  68 // executed.
  69 //----------------------------------------------------------------------
  70
  71 int
  72 l_fp_scmp(const l_fp first, const l_fp second)
  73 {
  74         u_int32 a[3], b[3];
  75
  76         const l_fp op1 = first;
  77         const l_fp op2 = second;
  78
  79         a[0] = op1.l_uf; a[1] = op1.l_ui; a[2] = 0;
  80         b[0] = op2.l_uf; b[1] = op2.l_ui; b[2] = 0;
  81
  82         a[2] -= (op1.l_i < 0);
  83         b[2] -= (op2.l_i < 0);
  84
  85         return cmp_work(a,b);
  86 }
  87
  88 int
  89 l_fp_ucmp(const l_fp first, l_fp second)
  90 {
  91         u_int32 a[3], b[3];
  92         const l_fp op1 = first;
  93         const l_fp op2 = second;
  94
  95         a[0] = op1.l_uf; a[1] = op1.l_ui; a[2] = 0;
  96         b[0] = op2.l_uf; b[1] = op2.l_ui; b[2] = 0;
  97
  98         return cmp_work(a,b);
  99 }
 100
 101 // maybe rename it to lf_cmp_work
 102 int
 103 cmp_work(u_int32 a[3], u_int32 b[3])
 104 {
 105         u_int32 cy, idx, tmp;
 106         for (cy = idx = 0; idx < 3; ++idx) {
 107                 tmp = a[idx]; cy  = (a[idx] -=   cy  ) > tmp;
 108                 tmp = a[idx]; cy |= (a[idx] -= b[idx]) > tmp;
 109         }
 110         if (a[2])
 111                 return -1;
 112         return a[0] || a[1];
 113 }
 114
 115
 116 //----------------------------------------------------------------------
 117 // imlementation of the LFP stuff
 118 // This should be easy enough...
 119 //----------------------------------------------------------------------
 120
 121 l_fp
 122 l_fp_init(int32 i, u_int32 f)
 123 {
 124         l_fp temp;
 125         temp.l_i  = i;
 126         temp.l_uf = f;
 127
 128         return temp;
 129 }
 130
 131 l_fp
 132 l_fp_add(const l_fp first, const l_fp second)
 133 {
 134         l_fp temp = first;
 135         L_ADD(&temp, &second);
 136
 137         return temp;
 138 }
 139
 140 l_fp
 141 l_fp_subtract(const l_fp first, const l_fp second)
 142 {
 143         l_fp temp = first;
 144         L_SUB(&temp, &second);
 145
 146         return temp;
 147 }
 148
 149 l_fp
 150 l_fp_negate(const l_fp first)
 151 {
 152         l_fp temp = first;
 153         L_NEG(&temp);
 154
 155         return temp;
 156 }
 157
 158 l_fp
 159 l_fp_abs(const l_fp first)
 160 {
 161         l_fp temp = first;
 162         if (L_ISNEG(&temp))
 163                 L_NEG(&temp);
 164         return temp;
 165 }
 166
 167 int
 168 l_fp_signum(const l_fp first)
 169 {
 170         if (first.l_ui & 0x80000000u)
 171                 return -1;
 172         return (first.l_ui || first.l_uf);
 173 }
 174
 175 double
 176 l_fp_convert_to_double(const l_fp first)
 177 {
 178         double res;
 179         LFPTOD(&first, res);
 180         return res;
 181 }
 182
 183 l_fp
 184 l_fp_init_from_double( double rhs)
 185 {
 186         l_fp temp;
 187         DTOLFP(rhs, &temp);
 188         return temp;
 189 }
 190
 191 void
 192 l_fp_swap(l_fp * first, l_fp *second)
 193 {
 194         l_fp temp = *second;
 195
 196         *second = *first;
 197         *first = temp;
 198
 199         return;
 200 }
 201
 202 //----------------------------------------------------------------------
 203 // testing the relational macros works better with proper predicate
 204 // formatting functions; it slows down the tests a bit, but makes for
 205 // readable failure messages.
 206 //----------------------------------------------------------------------
 207
 208
 209 bool
 210 l_isgt (const l_fp first, const l_fp second)
 211 {
 212
 213         return L_ISGT(&first, &second);
 214 }
 215
 216 bool
 217 l_isgtu(const l_fp first, const l_fp second)
 218 {
 219
 220         return L_ISGTU(&first, &second);
 221 }
 222
 223 bool
 224 l_ishis(const l_fp first, const l_fp second)
 225 {
 226
 227         return L_ISHIS(&first, &second);
 228 }
 229
 230 bool
 231 l_isgeq(const l_fp first, const l_fp second)
 232 {
 233
 234         return L_ISGEQ(&first, &second);
 235 }
 236
 237 bool
 238 l_isequ(const l_fp first, const l_fp second)
 239 {
 240
 241         return L_ISEQU(&first, &second);
 242 }
 243
 244
 245 //----------------------------------------------------------------------
 246 // test data table for add/sub and compare
 247 //----------------------------------------------------------------------
 248
 249
 250 static const lfp_hl addsub_tab[][3] = {
 251         // trivial idendity:
 252         {{0 ,0         }, { 0,0         }, { 0,0}},
 253         // with carry from fraction and sign change:
 254         {{-1,0x80000000}, { 0,0x80000000}, { 0,0}},
 255         // without carry from fraction
 256         {{ 1,0x40000000}, { 1,0x40000000}, { 2,0x80000000}},
 257         // with carry from fraction:
 258         {{ 1,0xC0000000}, { 1,0xC0000000}, { 3,0x80000000}},
 259         // with carry from fraction and sign change:
 260         {{0x7FFFFFFF, 0x7FFFFFFF}, {0x7FFFFFFF,0x7FFFFFFF}, {0xFFFFFFFE,0xFFFFFFFE}},
 261         // two tests w/o carry (used for l_fp<-->double):
 262         {{0x55555555,0xAAAAAAAA}, {0x11111111,0x11111111}, {0x66666666,0xBBBBBBBB}},
 263         {{0x55555555,0x55555555}, {0x11111111,0x11111111}, {0x66666666,0x66666666}},
 264         // wide-range test, triggers compare trouble
 265         {{0x80000000,0x00000001}, {0xFFFFFFFF,0xFFFFFFFE}, {0x7FFFFFFF,0xFFFFFFFF}}
 266 };
 267 static const size_t addsub_cnt = (sizeof(addsub_tab)/sizeof(addsub_tab[0]));
 268 static const size_t addsub_tot = (sizeof(addsub_tab)/sizeof(addsub_tab[0][0]));
 269
 270
 271
 272 //----------------------------------------------------------------------
 273 // epsilon estimation for the precision of a conversion double --> l_fp
 274 //
 275 // The error estimation limit is as follows:
 276 //  * The 'l_fp' fixed point fraction has 32 bits precision, so we allow
 277 //    for the LSB to toggle by clamping the epsilon to be at least 2^(-31)
 278 //
 279 //  * The double mantissa has a precsion 54 bits, so the other minimum is
 280 //    dval * (2^(-53))
 281 //
 282 //  The maximum of those two boundaries is used for the check.
 283 //
 284 // Note: once there are more than 54 bits between the highest and lowest
 285 // '1'-bit of the l_fp value, the roundtrip *will* create truncation
 286 // errors. This is an inherent property caused by the 54-bit mantissa of
 287 // the 'double' type.
 288 double
 289 eps(double d)
 290 {
 291
 292         return fmax(ldexp(1.0, -31), ldexp(fabs(d), -53));
 293 }
 294
 295 //----------------------------------------------------------------------
 296 // test addition
 297 //----------------------------------------------------------------------
 298 void
 299 test_AdditionLR(void)
 300 {
 301         size_t idx = 0;
 302
 303         for (idx = 0; idx < addsub_cnt; ++idx) {
 304                 l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 305                 l_fp op2 = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
 306                 l_fp e_res = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
 307                 l_fp res = l_fp_add(op1, op2);
 308
 309                 TEST_ASSERT_EQUAL_l_fp(e_res, res);
 310         }
 311         return;
 312 }
 313
 314 void
 315 test_AdditionRL(void)
 316 {
 317         size_t idx = 0;
 318
 319         for (idx = 0; idx < addsub_cnt; ++idx) {
 320                 l_fp op2 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 321                 l_fp op1 = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
 322                 l_fp e_res = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
 323                 l_fp res = l_fp_add(op1, op2);
 324
 325                 TEST_ASSERT_EQUAL_l_fp(e_res, res);
 326         }
 327         return;
 328 }
 329
 330
 331 //----------------------------------------------------------------------
 332 // test subtraction
 333 //----------------------------------------------------------------------
 334 void
 335 test_SubtractionLR(void)
 336 {
 337         size_t idx = 0;
 338
 339         for (idx = 0; idx < addsub_cnt; ++idx) {
 340                 l_fp op2 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 341                 l_fp e_res = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
 342                 l_fp op1 = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
 343                 l_fp res = l_fp_subtract(op1, op2);
 344
 345                 TEST_ASSERT_EQUAL_l_fp(e_res, res);
 346         }
 347         return;
 348 }
 349
 350 void
 351 test_SubtractionRL(void)
 352 {
 353         size_t idx = 0;
 354
 355         for (idx = 0; idx < addsub_cnt; ++idx) {
 356                 l_fp e_res = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 357                 l_fp op2 = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
 358                 l_fp op1 = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
 359                 l_fp res = l_fp_subtract(op1, op2);
 360
 361                 TEST_ASSERT_EQUAL_l_fp(e_res, res);
 362         }
 363         return;
 364 }
 365
 366 //----------------------------------------------------------------------
 367 // test negation
 368 //----------------------------------------------------------------------
 369
 370 void
 371 test_Negation(void)
 372 {
 373         size_t idx = 0;
 374
 375         for (idx = 0; idx < addsub_cnt; ++idx) {
 376                 l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 377                 l_fp op2 = l_fp_negate(op1);
 378                 l_fp sum = l_fp_add(op1, op2);
 379
 380                 l_fp zero = l_fp_init(0, 0);
 381
 382                 TEST_ASSERT_EQUAL_l_fp(zero, sum);
 383         }
 384         return;
 385 }
 386
 387
 388
 389 //----------------------------------------------------------------------
 390 // test absolute value
 391 //----------------------------------------------------------------------
 392 void
 393 test_Absolute(void)
 394 {
 395         size_t idx = 0;
 396
 397         for (idx = 0; idx < addsub_cnt; ++idx) {
 398                 l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 399                 l_fp op2 = l_fp_abs(op1);
 400
 401                 TEST_ASSERT_TRUE(l_fp_signum(op2) >= 0);
 402
 403                 if (l_fp_signum(op1) >= 0)
 404                         op1 = l_fp_subtract(op1, op2);
 405                 else
 406                         op1 = l_fp_add(op1, op2);
 407
 408                 l_fp zero = l_fp_init(0, 0);
 409
 410                 TEST_ASSERT_EQUAL_l_fp(zero, op1);
 411         }
 412
 413         // There is one special case we have to check: the minimum
 414         // value cannot be negated, or, to be more precise, the
 415         // negation reproduces the original pattern.
 416         l_fp minVal = l_fp_init(0x80000000, 0x00000000);
 417         l_fp minAbs = l_fp_abs(minVal);
 418         TEST_ASSERT_EQUAL(-1, l_fp_signum(minVal));
 419
 420         TEST_ASSERT_EQUAL_l_fp(minVal, minAbs);
 421
 422         return;
 423 }
 424
 425
 426 //----------------------------------------------------------------------
 427 // fp -> double -> fp rountrip test
 428 //----------------------------------------------------------------------
 429 void
 430 test_FDF_RoundTrip(void)
 431 {
 432         size_t idx = 0;
 433
 434         // since a l_fp has 64 bits in it's mantissa and a double has
 435         // only 54 bits available (including the hidden '1') we have to
 436         // make a few concessions on the roundtrip precision. The 'eps()'
 437         // function makes an educated guess about the avilable precision
 438         // and checks the difference in the two 'l_fp' values against
 439         // that limit.
 440
 441         for (idx = 0; idx < addsub_cnt; ++idx) {
 442                 l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 443                 double op2 = l_fp_convert_to_double(op1);
 444                 l_fp op3 = l_fp_init_from_double(op2);
 445
 446                 l_fp temp = l_fp_subtract(op1, op3);
 447                 double d = l_fp_convert_to_double(temp);
 448                 TEST_ASSERT_DOUBLE_WITHIN(eps(op2), 0.0, fabs(d));
 449         }
 450
 451         return;
 452 }
 453
 454
 455 //----------------------------------------------------------------------
 456 // test the compare stuff
 457 //
 458 // This uses the local compare and checks if the operations using the
 459 // macros in 'ntp_fp.h' produce mathing results.
 460 // ----------------------------------------------------------------------
 461 void
 462 test_SignedRelOps(void)
 463 {
 464         const lfp_hl * tv = (&addsub_tab[0][0]);
 465         size_t lc ;
 466
 467         for (lc = addsub_tot - 1; lc; --lc, ++tv) {
 468                 l_fp op1 = l_fp_init(tv[0].h, tv[0].l);
 469                 l_fp op2 = l_fp_init(tv[1].h, tv[1].l);
 470                 int cmp = l_fp_scmp(op1, op2);
 471
 472                 switch (cmp) {
 473                 case -1:
 474                         //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
 475                         l_fp_swap(&op1, &op2);
 476                         //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
 477                 case 1:
 478                         TEST_ASSERT_TRUE (l_isgt(op1, op2));
 479                         TEST_ASSERT_FALSE(l_isgt(op2, op1));
 480
 481                         TEST_ASSERT_TRUE (l_isgeq(op1, op2));
 482                         TEST_ASSERT_FALSE(l_isgeq(op2, op1));
 483
 484                         TEST_ASSERT_FALSE(l_isequ(op1, op2));
 485                         TEST_ASSERT_FALSE(l_isequ(op2, op1));
 486                         break;
 487                 case 0:
 488                         TEST_ASSERT_FALSE(l_isgt(op1, op2));
 489                         TEST_ASSERT_FALSE(l_isgt(op2, op1));
 490
 491                         TEST_ASSERT_TRUE (l_isgeq(op1, op2));
 492                         TEST_ASSERT_TRUE (l_isgeq(op2, op1));
 493
 494                         TEST_ASSERT_TRUE (l_isequ(op1, op2));
 495                         TEST_ASSERT_TRUE (l_isequ(op2, op1));
 496                         break;
 497                 default:
 498                         TEST_FAIL_MESSAGE("unexpected UCMP result: ");
 499                 }
 500         }
 501
 502         return;
 503 }
 504
 505 void
 506 test_UnsignedRelOps(void)
 507 {
 508         const lfp_hl * tv =(&addsub_tab[0][0]);
 509         size_t lc;
 510
 511         for (lc = addsub_tot - 1; lc; --lc, ++tv) {
 512                 l_fp op1 = l_fp_init(tv[0].h, tv[0].l);
 513                 l_fp op2 = l_fp_init(tv[1].h, tv[1].l);
 514                 int cmp = l_fp_ucmp(op1, op2);
 515
 516                 switch (cmp) {
 517                 case -1:
 518                         //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
 519                         l_fp_swap(&op1, &op2);
 520                         //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
 521                 case 1:
 522                         TEST_ASSERT_TRUE (l_isgtu(op1, op2));
 523                         TEST_ASSERT_FALSE(l_isgtu(op2, op1));
 524
 525                         TEST_ASSERT_TRUE (l_ishis(op1, op2));
 526                         TEST_ASSERT_FALSE(l_ishis(op2, op1));
 527                         break;
 528                 case 0:
 529                         TEST_ASSERT_FALSE(l_isgtu(op1, op2));
 530                         TEST_ASSERT_FALSE(l_isgtu(op2, op1));
 531
 532                         TEST_ASSERT_TRUE (l_ishis(op1, op2));
 533                         TEST_ASSERT_TRUE (l_ishis(op2, op1));
 534                         break;
 535                 default:
 536                         TEST_FAIL_MESSAGE("unexpected UCMP result: ");
 537                 }
 538         }
 539
 540         return;
 541 }
 542
 543 /*
 544 */
 545
 546 //----------------------------------------------------------------------
 547 // that's all folks... but feel free to add things!
 548 //----------------------------------------------------------------------