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 /* replaced TEST_ASSERT_EQUAL_MEMORY(&a, &b, sizeof(a)) with TEST_ASSERT_EQUAL_l_fp(a, b).
  13    It's safer this way, because structs can be compared even if they aren't initiated
  14    with memset (due to padding bytes).
  15 */
  16 #define TEST_ASSERT_EQUAL_l_fp(a, b) { \
  17     TEST_ASSERT_EQUAL_MESSAGE(a.l_i, b.l_i, "Field l_i"); \
  18     TEST_ASSERT_EQUAL_UINT_MESSAGE(a.l_uf, b.l_uf, "Field l_uf");       \
  19 }
  20
  21
  22
  23 typedef int bool; // typedef enum { FALSE, TRUE } boolean; -> can't use this because TRUE and FALSE are already defined
  24
  25
  26 typedef struct  {
  27         uint32_t h, l;
  28 } lfp_hl;
  29
  30
  31 int l_fp_scmp(const l_fp first, const l_fp second);
  32 int l_fp_ucmp(const l_fp first, l_fp second );
  33 l_fp l_fp_init(int32 i, u_int32 f);
  34 l_fp l_fp_add(const l_fp first, const l_fp second);
  35 l_fp l_fp_subtract(const l_fp first, const l_fp second);
  36 l_fp l_fp_negate(const l_fp first);
  37 l_fp l_fp_abs(const l_fp first);
  38 int l_fp_signum(const l_fp first);
  39 double l_fp_convert_to_double(const l_fp first);
  40 l_fp l_fp_init_from_double( double rhs);
  41 void l_fp_swap(l_fp * first, l_fp *second);
  42 bool l_isgt(const l_fp first, const l_fp second);
  43 bool l_isgtu(const l_fp first, const l_fp second);
  44 bool l_ishis(const l_fp first, const l_fp second);
  45 bool l_isgeq(const l_fp first, const l_fp second);
  46 bool l_isequ(const l_fp first, const l_fp second);
  47 double eps(double d);
  48
  49
  50 void test_AdditionLR(void);
  51 void test_AdditionRL(void);
  52 void test_SubtractionLR(void);
  53 void test_SubtractionRL(void);
  54 void test_Negation(void);
  55 void test_Absolute(void);
  56 void test_FDF_RoundTrip(void);
  57 void test_SignedRelOps(void);
  58 void test_UnsignedRelOps(void);
  59
  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         l_fp temp = *second;
 194
 195         *second = *first;
 196         *first = temp;
 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         return L_ISGT(&first, &second);
 209 }
 210
 211 bool
 212 l_isgtu(const l_fp first, const l_fp second) {
 213         return L_ISGTU(&first, &second);
 214 }
 215
 216 bool
 217 l_ishis(const l_fp first, const l_fp second) {
 218         return L_ISHIS(&first, &second);
 219 }
 220
 221 bool
 222 l_isgeq(const l_fp first, const l_fp second) {
 223         return L_ISGEQ(&first, &second);
 224 }
 225
 226 bool
 227 l_isequ(const l_fp first, const l_fp second) {
 228         return L_ISEQU(&first, &second);
 229 }
 230
 231
 232 //----------------------------------------------------------------------
 233 // test data table for add/sub and compare
 234 //----------------------------------------------------------------------
 235
 236
 237 static const lfp_hl addsub_tab[][3] = {
 238         // trivial idendity:
 239         {{0 ,0         }, { 0,0         }, { 0,0}},
 240         // with carry from fraction and sign change:
 241         {{-1,0x80000000}, { 0,0x80000000}, { 0,0}},
 242         // without carry from fraction
 243         {{ 1,0x40000000}, { 1,0x40000000}, { 2,0x80000000}},
 244         // with carry from fraction:
 245         {{ 1,0xC0000000}, { 1,0xC0000000}, { 3,0x80000000}},
 246         // with carry from fraction and sign change:
 247         {{0x7FFFFFFF, 0x7FFFFFFF}, {0x7FFFFFFF,0x7FFFFFFF}, {0xFFFFFFFE,0xFFFFFFFE}},
 248         // two tests w/o carry (used for l_fp<-->double):
 249         {{0x55555555,0xAAAAAAAA}, {0x11111111,0x11111111}, {0x66666666,0xBBBBBBBB}},
 250         {{0x55555555,0x55555555}, {0x11111111,0x11111111}, {0x66666666,0x66666666}},
 251         // wide-range test, triggers compare trouble
 252         {{0x80000000,0x00000001}, {0xFFFFFFFF,0xFFFFFFFE}, {0x7FFFFFFF,0xFFFFFFFF}}
 253 };
 254 static const size_t addsub_cnt = (sizeof(addsub_tab)/sizeof(addsub_tab[0]));
 255 static const size_t addsub_tot = (sizeof(addsub_tab)/sizeof(addsub_tab[0][0]));
 256
 257
 258
 259 //----------------------------------------------------------------------
 260 // epsilon estimation for the precision of a conversion double --> l_fp
 261 //
 262 // The error estimation limit is as follows:
 263 //  * The 'l_fp' fixed point fraction has 32 bits precision, so we allow
 264 //    for the LSB to toggle by clamping the epsilon to be at least 2^(-31)
 265 //
 266 //  * The double mantissa has a precsion 54 bits, so the other minimum is
 267 //    dval * (2^(-53))
 268 //
 269 //  The maximum of those two boundaries is used for the check.
 270 //
 271 // Note: once there are more than 54 bits between the highest and lowest
 272 // '1'-bit of the l_fp value, the roundtrip *will* create truncation
 273 // errors. This is an inherent property caused by the 54-bit mantissa of
 274 // the 'double' type.
 275 double
 276 eps(double d)
 277 {
 278         return fmax(ldexp(1.0, -31), ldexp(fabs(d), -53));
 279 }
 280
 281 //----------------------------------------------------------------------
 282 // test addition
 283 //----------------------------------------------------------------------
 284 void
 285 test_AdditionLR(void) {
 286
 287         size_t idx = 0;
 288         for (idx = 0; idx < addsub_cnt; ++idx) {
 289                 l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 290                 l_fp op2 = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
 291                 l_fp exp = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
 292                 l_fp res = l_fp_add(op1, op2);
 293
 294                 TEST_ASSERT_EQUAL_l_fp(exp, res);
 295         }
 296 }
 297
 298 void
 299 test_AdditionRL(void) {
 300         size_t idx = 0;
 301         for (idx = 0; idx < addsub_cnt; ++idx) {
 302                 l_fp op2 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 303                 l_fp op1 = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
 304                 l_fp exp = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
 305                 l_fp res = l_fp_add(op1, op2);
 306
 307                 TEST_ASSERT_EQUAL_l_fp(exp, res);
 308         }
 309 }
 310
 311
 312
 313 //----------------------------------------------------------------------
 314 // test subtraction
 315 //----------------------------------------------------------------------
 316 void
 317 test_SubtractionLR(void) {
 318         size_t idx = 0;
 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 exp = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
 322                 l_fp op1 = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
 323                 l_fp res = l_fp_subtract(op1, op2);
 324
 325                 TEST_ASSERT_EQUAL_l_fp(exp, res);
 326         }
 327 }
 328
 329 void
 330 test_SubtractionRL(void) {
 331         size_t idx = 0;
 332         for (idx = 0; idx < addsub_cnt; ++idx) {
 333                 l_fp exp = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 334                 l_fp op2 = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
 335                 l_fp op1 = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
 336                 l_fp res = l_fp_subtract(op1, op2);
 337
 338                 TEST_ASSERT_EQUAL_l_fp(exp, res);
 339         }
 340 }
 341
 342 //----------------------------------------------------------------------
 343 // test negation
 344 //----------------------------------------------------------------------
 345
 346 void
 347 test_Negation(void) {
 348
 349         size_t idx = 0;
 350         for (idx = 0; idx < addsub_cnt; ++idx) {
 351                 l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 352                 l_fp op2 = l_fp_negate(op1);
 353                 l_fp sum = l_fp_add(op1, op2);
 354
 355                 l_fp zero = l_fp_init(0, 0);
 356
 357                 TEST_ASSERT_EQUAL_l_fp(zero, sum);
 358         }
 359 }
 360
 361
 362
 363 //----------------------------------------------------------------------
 364 // test absolute value
 365 //----------------------------------------------------------------------
 366 void
 367 test_Absolute(void) {
 368         size_t idx = 0;
 369         for (idx = 0; idx < addsub_cnt; ++idx) {
 370                 l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 371                 l_fp op2 = l_fp_abs(op1);
 372
 373                 TEST_ASSERT_TRUE(l_fp_signum(op2) >= 0);
 374
 375                 if (l_fp_signum(op1) >= 0)
 376                         op1 = l_fp_subtract(op1, op2);
 377                 else
 378                         op1 = l_fp_add(op1, op2);
 379
 380                 l_fp zero = l_fp_init(0, 0);
 381
 382                 TEST_ASSERT_EQUAL_l_fp(zero, op1);
 383         }
 384
 385         // There is one special case we have to check: the minimum
 386         // value cannot be negated, or, to be more precise, the
 387         // negation reproduces the original pattern.
 388         l_fp minVal = l_fp_init(0x80000000, 0x00000000);
 389         l_fp minAbs = l_fp_abs(minVal);
 390         TEST_ASSERT_EQUAL(-1, l_fp_signum(minVal));
 391
 392         TEST_ASSERT_EQUAL_l_fp(minVal, minAbs);
 393 }
 394
 395
 396 //----------------------------------------------------------------------
 397 // fp -> double -> fp rountrip test
 398 //----------------------------------------------------------------------
 399 void
 400 test_FDF_RoundTrip(void) {
 401         // since a l_fp has 64 bits in it's mantissa and a double has
 402         // only 54 bits available (including the hidden '1') we have to
 403         // make a few concessions on the roundtrip precision. The 'eps()'
 404         // function makes an educated guess about the avilable precision
 405         // and checks the difference in the two 'l_fp' values against
 406         // that limit.
 407         size_t idx = 0;
 408         for (idx = 0; idx < addsub_cnt; ++idx) {
 409                 l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
 410                 double op2 = l_fp_convert_to_double(op1);
 411                 l_fp op3 = l_fp_init_from_double(op2);
 412
 413                 l_fp temp = l_fp_subtract(op1, op3);
 414                 double d = l_fp_convert_to_double(temp);
 415                 TEST_ASSERT_DOUBLE_WITHIN(eps(op2), 0.0, fabs(d));
 416         }
 417 }
 418
 419
 420 //----------------------------------------------------------------------
 421 // test the compare stuff
 422 //
 423 // This uses the local compare and checks if the operations using the
 424 // macros in 'ntp_fp.h' produce mathing results.
 425 // ----------------------------------------------------------------------
 426 void
 427 test_SignedRelOps(void) {
 428         const lfp_hl * tv = (&addsub_tab[0][0]);
 429         size_t lc ;
 430         for (lc = addsub_tot - 1; lc; --lc, ++tv) {
 431                 l_fp op1 = l_fp_init(tv[0].h, tv[0].l);
 432                 l_fp op2 = l_fp_init(tv[1].h, tv[1].l);
 433                 int cmp = l_fp_scmp(op1, op2);
 434
 435                 switch (cmp) {
 436                 case -1:
 437                         //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
 438                         l_fp_swap(&op1, &op2);
 439                         //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
 440                 case 1:
 441                         TEST_ASSERT_TRUE (l_isgt(op1, op2));
 442                         TEST_ASSERT_FALSE(l_isgt(op2, op1));
 443
 444                         TEST_ASSERT_TRUE (l_isgeq(op1, op2));
 445                         TEST_ASSERT_FALSE(l_isgeq(op2, op1));
 446
 447                         TEST_ASSERT_FALSE(l_isequ(op1, op2));
 448                         TEST_ASSERT_FALSE(l_isequ(op2, op1));
 449                         break;
 450                 case 0:
 451                         TEST_ASSERT_FALSE(l_isgt(op1, op2));
 452                         TEST_ASSERT_FALSE(l_isgt(op2, op1));
 453
 454                         TEST_ASSERT_TRUE (l_isgeq(op1, op2));
 455                         TEST_ASSERT_TRUE (l_isgeq(op2, op1));
 456
 457                         TEST_ASSERT_TRUE (l_isequ(op1, op2));
 458                         TEST_ASSERT_TRUE (l_isequ(op2, op1));
 459                         break;
 460                 default:
 461                         TEST_FAIL_MESSAGE("unexpected UCMP result: " );
 462                 }
 463         }
 464 }
 465
 466 void
 467 test_UnsignedRelOps(void) {
 468         const lfp_hl * tv =(&addsub_tab[0][0]);
 469         size_t lc;
 470         for (lc = addsub_tot - 1; lc; --lc, ++tv) {
 471                 l_fp op1 = l_fp_init(tv[0].h, tv[0].l);
 472                 l_fp op2 = l_fp_init(tv[1].h, tv[1].l);
 473                 int cmp = l_fp_ucmp(op1, op2);
 474
 475                 switch (cmp) {
 476                 case -1:
 477                         //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
 478                         l_fp_swap(&op1, &op2);
 479                         //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
 480                 case 1:
 481                         TEST_ASSERT_TRUE (l_isgtu(op1, op2));
 482                         TEST_ASSERT_FALSE(l_isgtu(op2, op1));
 483
 484                         TEST_ASSERT_TRUE (l_ishis(op1, op2));
 485                         TEST_ASSERT_FALSE(l_ishis(op2, op1));
 486                         break;
 487                 case 0:
 488                         TEST_ASSERT_FALSE(l_isgtu(op1, op2));
 489                         TEST_ASSERT_FALSE(l_isgtu(op2, op1));
 490
 491                         TEST_ASSERT_TRUE (l_ishis(op1, op2));
 492                         TEST_ASSERT_TRUE (l_ishis(op2, op1));
 493                         break;
 494                 default:
 495                         TEST_FAIL_MESSAGE("unexpected UCMP result: " );
 496                 }
 497         }
 498 }
 499 /*
 500 */
 501
 502 //----------------------------------------------------------------------
 503 // that's all folks... but feel free to add things!
 504 //----------------------------------------------------------------------