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 26typedef struct { 27 uint32_t h, l; 28} lfp_hl; 29 30 31int l_fp_scmp(const l_fp first, const l_fp second); 32int l_fp_ucmp(const l_fp first, l_fp second); 33l_fp l_fp_init(int32 i, u_int32 f); 34l_fp l_fp_add(const l_fp first, const l_fp second); 35l_fp l_fp_subtract(const l_fp first, const l_fp second); 36l_fp l_fp_negate(const l_fp first); 37l_fp l_fp_abs(const l_fp first); 38int l_fp_signum(const l_fp first); 39double l_fp_convert_to_double(const l_fp first); 40l_fp l_fp_init_from_double( double rhs); 41void l_fp_swap(l_fp * first, l_fp *second); 42bool l_isgt(const l_fp first, const l_fp second); 43bool l_isgtu(const l_fp first, const l_fp second); 44bool l_ishis(const l_fp first, const l_fp second); 45bool l_isgeq(const l_fp first, const l_fp second); 46bool l_isequ(const l_fp first, const l_fp second); 47double eps(double d); 48 49 50void test_AdditionLR(void); 51void test_AdditionRL(void); 52void test_SubtractionLR(void); 53void test_SubtractionRL(void); 54void test_Negation(void); 55void test_Absolute(void); 56void test_FDF_RoundTrip(void); 57void test_SignedRelOps(void); 58void test_UnsignedRelOps(void); 59 60 61static 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 70int 71l_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 87int 88l_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 101int 102cmp_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 120l_fp 121l_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 130l_fp 131l_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 139l_fp 140l_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 148l_fp 149l_fp_negate(const l_fp first) 150{ 151 l_fp temp = first; 152 L_NEG(&temp); 153 154 return temp; 155} 156 157l_fp 158l_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 166int 167l_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 174double 175l_fp_convert_to_double(const l_fp first) 176{ 177 double res; 178 LFPTOD(&first, res); 179 return res; 180} 181 182l_fp 183l_fp_init_from_double( double rhs) 184{ 185 l_fp temp; 186 DTOLFP(rhs, &temp); 187 return temp; 188} 189 190void 191l_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 208bool 209l_isgt (const l_fp first, const l_fp second) 210{ 211 212 return L_ISGT(&first, &second); 213} 214 215bool 216l_isgtu(const l_fp first, const l_fp second) 217{ 218 219 return L_ISGTU(&first, &second); 220} 221 222bool 223l_ishis(const l_fp first, const l_fp second) 224{ 225 226 return L_ISHIS(&first, &second); 227} 228 229bool 230l_isgeq(const l_fp first, const l_fp second) 231{ 232 233 return L_ISGEQ(&first, &second); 234} 235 236bool 237l_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 249static 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}; 266static const size_t addsub_cnt = (sizeof(addsub_tab)/sizeof(addsub_tab[0])); 267static 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. 287double 288eps(double d) 289{ 290 291 return fmax(ldexp(1.0, -31), ldexp(fabs(d), -53)); 292} 293 294//---------------------------------------------------------------------- 295// test addition 296//---------------------------------------------------------------------- 297void 298test_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 313void 314test_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//---------------------------------------------------------------------- 333void 334test_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 349void 350test_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 369void 370test_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//---------------------------------------------------------------------- 391void 392test_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//---------------------------------------------------------------------- 428void 429test_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// ---------------------------------------------------------------------- 460void 461test_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 504void 505test_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