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 //---------------------------------------------------------------------- 546