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