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
l_fp_scmp(const l_fp first,const l_fp second)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
l_fp_ucmp(const l_fp first,l_fp second)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
cmp_work(u_int32 a[3],u_int32 b[3])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
l_fp_init(int32 i,u_int32 f)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
l_fp_add(const l_fp first,const l_fp second)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
l_fp_subtract(const l_fp first,const l_fp second)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
l_fp_negate(const l_fp first)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
l_fp_abs(const l_fp first)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
l_fp_signum(const l_fp first)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
l_fp_convert_to_double(const l_fp first)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
l_fp_init_from_double(double rhs)183 l_fp_init_from_double( double rhs)
184 {
185 l_fp temp;
186 DTOLFP(rhs, &temp);
187 return temp;
188 }
189
190 void
l_fp_swap(l_fp * first,l_fp * second)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
l_isgt(const l_fp first,const l_fp second)209 l_isgt (const l_fp first, const l_fp second)
210 {
211
212 return L_ISGT(&first, &second);
213 }
214
215 bool
l_isgtu(const l_fp first,const l_fp second)216 l_isgtu(const l_fp first, const l_fp second)
217 {
218
219 return L_ISGTU(&first, &second);
220 }
221
222 bool
l_ishis(const l_fp first,const l_fp second)223 l_ishis(const l_fp first, const l_fp second)
224 {
225
226 return L_ISHIS(&first, &second);
227 }
228
229 bool
l_isgeq(const l_fp first,const l_fp second)230 l_isgeq(const l_fp first, const l_fp second)
231 {
232
233 return L_ISGEQ(&first, &second);
234 }
235
236 bool
l_isequ(const l_fp first,const l_fp second)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
eps(double d)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
test_AdditionLR(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
test_AdditionRL(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
test_SubtractionLR(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
test_SubtractionRL(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
test_Negation(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
test_Absolute(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
test_FDF_RoundTrip(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
test_SignedRelOps(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
test_UnsignedRelOps(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