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