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