lfpfunc.c revision 284990
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 19typedef struct { 20 uint32_t h, l; 21} lfp_hl; 22 23 24static 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 71static 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//---------------------------------------------------------------------- 82int 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 99int 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 ??? 113int 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 133l_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 144l_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 152l_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 161l_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 170l_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 178int 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 185double l_fp_convert_to_double(const l_fp first) 186{ 187 double res; 188 LFPTOD(&first, res); 189 return res; 190} 191 192l_fp l_fp_init_from_double( double rhs) 193{ 194 l_fp temp; 195 DTOLFP(rhs, &temp); 196 return temp; 197} 198 199 200 201void 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/* 211LFP::LFP() 212{ 213 _v.l_ui = 0; 214 _v.l_uf = 0; 215} 216 217 218 219std::string 220LFP::toString() const 221{ 222 std::ostringstream oss; 223 toStream(oss); 224 return oss.str(); 225} 226 227std::ostream& 228LFP::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 237bool 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 253typedef int bool; //typedef enum { FALSE, TRUE } boolean; -> can't use this because TRUE and FALSE are already defined 254 255 256bool l_isgt (const l_fp first, const l_fp second) 257 { return L_ISGT(&first, &second); } 258bool l_isgtu(const l_fp first, const l_fp second) 259 { return L_ISGTU(&first, &second); } 260bool l_ishis(const l_fp first, const l_fp second) 261 { return L_ISHIS(&first, &second); } 262bool l_isgeq(const l_fp first, const l_fp second) 263 { return L_ISGEQ(&first, &second); } 264bool 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 273static 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}; 290static const size_t addsub_cnt = (sizeof(addsub_tab)/sizeof(addsub_tab[0])); 291static 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. 311double 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//---------------------------------------------------------------------- 321void 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 340void 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//---------------------------------------------------------------------- 359void 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 374void 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 392void 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//---------------------------------------------------------------------- 413void 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//---------------------------------------------------------------------- 448void 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// ---------------------------------------------------------------------- 480void 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 523void 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//---------------------------------------------------------------------- 562