lfpfunc.c revision 293896
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 24typedef int bool; // typedef enum { FALSE, TRUE } boolean; -> can't use this because TRUE and FALSE are already defined 25 26 27typedef struct { 28 uint32_t h, l; 29} lfp_hl; 30 31 32int l_fp_scmp(const l_fp first, const l_fp second); 33int l_fp_ucmp(const l_fp first, l_fp second); 34l_fp l_fp_init(int32 i, u_int32 f); 35l_fp l_fp_add(const l_fp first, const l_fp second); 36l_fp l_fp_subtract(const l_fp first, const l_fp second); 37l_fp l_fp_negate(const l_fp first); 38l_fp l_fp_abs(const l_fp first); 39int l_fp_signum(const l_fp first); 40double l_fp_convert_to_double(const l_fp first); 41l_fp l_fp_init_from_double( double rhs); 42void l_fp_swap(l_fp * first, l_fp *second); 43bool l_isgt(const l_fp first, const l_fp second); 44bool l_isgtu(const l_fp first, const l_fp second); 45bool l_ishis(const l_fp first, const l_fp second); 46bool l_isgeq(const l_fp first, const l_fp second); 47bool l_isequ(const l_fp first, const l_fp second); 48double eps(double d); 49 50 51void test_AdditionLR(void); 52void test_AdditionRL(void); 53void test_SubtractionLR(void); 54void test_SubtractionRL(void); 55void test_Negation(void); 56void test_Absolute(void); 57void test_FDF_RoundTrip(void); 58void test_SignedRelOps(void); 59void test_UnsignedRelOps(void); 60 61 62static 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 71int 72l_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 88int 89l_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 102int 103cmp_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 121l_fp 122l_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 131l_fp 132l_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 140l_fp 141l_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 149l_fp 150l_fp_negate(const l_fp first) 151{ 152 l_fp temp = first; 153 L_NEG(&temp); 154 155 return temp; 156} 157 158l_fp 159l_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 167int 168l_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 175double 176l_fp_convert_to_double(const l_fp first) 177{ 178 double res; 179 LFPTOD(&first, res); 180 return res; 181} 182 183l_fp 184l_fp_init_from_double( double rhs) 185{ 186 l_fp temp; 187 DTOLFP(rhs, &temp); 188 return temp; 189} 190 191void 192l_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 209bool 210l_isgt (const l_fp first, const l_fp second) 211{ 212 213 return L_ISGT(&first, &second); 214} 215 216bool 217l_isgtu(const l_fp first, const l_fp second) 218{ 219 220 return L_ISGTU(&first, &second); 221} 222 223bool 224l_ishis(const l_fp first, const l_fp second) 225{ 226 227 return L_ISHIS(&first, &second); 228} 229 230bool 231l_isgeq(const l_fp first, const l_fp second) 232{ 233 234 return L_ISGEQ(&first, &second); 235} 236 237bool 238l_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 250static 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}; 267static const size_t addsub_cnt = (sizeof(addsub_tab)/sizeof(addsub_tab[0])); 268static 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. 288double 289eps(double d) 290{ 291 292 return fmax(ldexp(1.0, -31), ldexp(fabs(d), -53)); 293} 294 295//---------------------------------------------------------------------- 296// test addition 297//---------------------------------------------------------------------- 298void 299test_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 314void 315test_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//---------------------------------------------------------------------- 334void 335test_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 350void 351test_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 370void 371test_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//---------------------------------------------------------------------- 392void 393test_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//---------------------------------------------------------------------- 429void 430test_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// ---------------------------------------------------------------------- 461void 462test_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 505void 506test_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//---------------------------------------------------------------------- 549