contrib/ntp/tests/libntp/lfpfunc.c

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