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