/*      $NetBSD: lfpfunc.c,v 1.2 2020/05/25 20:47:36 christos Exp $     */

#include "config.h"

#include "ntp_stdlib.h"
#include "ntp_fp.h"

#include "unity.h"

#include <float.h>
#include <math.h>


/*
  replaced:    TEST_ASSERT_EQUAL_MEMORY(&a, &b, sizeof(a))
  with:        TEST_ASSERT_EQUAL_l_fp(a, b).
  It's safer this way, because structs can be compared even if they
  aren't initiated with memset (due to padding bytes).
*/
#define TEST_ASSERT_EQUAL_l_fp(a, b) {                                  \
       TEST_ASSERT_EQUAL_MESSAGE(a.l_i, b.l_i, "Field l_i");           \
       TEST_ASSERT_EQUAL_UINT_MESSAGE(a.l_uf, b.l_uf, "Field l_uf");   \
}


typedef struct  {
       uint32_t h, l;
} lfp_hl;


int     l_fp_scmp(const l_fp first, const l_fp second);
int     l_fp_ucmp(const l_fp first, l_fp second);
l_fp    l_fp_init(int32 i, u_int32 f);
l_fp    l_fp_add(const l_fp first, const l_fp second);
l_fp    l_fp_subtract(const l_fp first, const l_fp second);
l_fp    l_fp_negate(const l_fp first);
l_fp    l_fp_abs(const l_fp first);
int     l_fp_signum(const l_fp first);
double  l_fp_convert_to_double(const l_fp first);
l_fp    l_fp_init_from_double( double rhs);
void    l_fp_swap(l_fp * first, l_fp *second);
bool    l_isgt(const l_fp first, const l_fp second);
bool    l_isgtu(const l_fp first, const l_fp second);
bool    l_ishis(const l_fp first, const l_fp second);
bool    l_isgeq(const l_fp first, const l_fp second);
bool    l_isequ(const l_fp first, const l_fp second);
double  eps(double d);


void test_AdditionLR(void);
void test_AdditionRL(void);
void test_SubtractionLR(void);
void test_SubtractionRL(void);
void test_Negation(void);
void test_Absolute(void);
void test_FDF_RoundTrip(void);
void test_SignedRelOps(void);
void test_UnsignedRelOps(void);


static int cmp_work(u_int32 a[3], u_int32 b[3]);

//----------------------------------------------------------------------
// reference comparision
// This is implementad as a full signed MP-subtract in 3 limbs, where
// the operands are zero or sign extended before the subtraction is
// executed.
//----------------------------------------------------------------------

int
l_fp_scmp(const l_fp first, const l_fp second)
{
       u_int32 a[3], b[3];

       const l_fp op1 = first;
       const l_fp op2 = second;

       a[0] = op1.l_uf; a[1] = op1.l_ui; a[2] = 0;
       b[0] = op2.l_uf; b[1] = op2.l_ui; b[2] = 0;

       a[2] -= (op1.l_i < 0);
       b[2] -= (op2.l_i < 0);

       return cmp_work(a,b);
}

int
l_fp_ucmp(const l_fp first, l_fp second)
{
       u_int32 a[3], b[3];
       const l_fp op1 = first;
       const l_fp op2 = second;

       a[0] = op1.l_uf; a[1] = op1.l_ui; a[2] = 0;
       b[0] = op2.l_uf; b[1] = op2.l_ui; b[2] = 0;

       return cmp_work(a,b);
}

// maybe rename it to lf_cmp_work
int
cmp_work(u_int32 a[3], u_int32 b[3])
{
       u_int32 cy, idx, tmp;
       for (cy = idx = 0; idx < 3; ++idx) {
               tmp = a[idx]; cy  = (a[idx] -=   cy  ) > tmp;
               tmp = a[idx]; cy |= (a[idx] -= b[idx]) > tmp;
       }
       if (a[2])
               return -1;
       return a[0] || a[1];
}


//----------------------------------------------------------------------
// imlementation of the LFP stuff
// This should be easy enough...
//----------------------------------------------------------------------

l_fp
l_fp_init(int32 i, u_int32 f)
{
       l_fp temp;
       temp.l_i  = i;
       temp.l_uf = f;

       return temp;
}

l_fp
l_fp_add(const l_fp first, const l_fp second)
{
       l_fp temp = first;
       L_ADD(&temp, &second);

       return temp;
}

l_fp
l_fp_subtract(const l_fp first, const l_fp second)
{
       l_fp temp = first;
       L_SUB(&temp, &second);

       return temp;
}

l_fp
l_fp_negate(const l_fp first)
{
       l_fp temp = first;
       L_NEG(&temp);

       return temp;
}

l_fp
l_fp_abs(const l_fp first)
{
       l_fp temp = first;
       if (L_ISNEG(&temp))
               L_NEG(&temp);
       return temp;
}

int
l_fp_signum(const l_fp first)
{
       if (first.l_ui & 0x80000000u)
               return -1;
       return (first.l_ui || first.l_uf);
}

double
l_fp_convert_to_double(const l_fp first)
{
       double res;
       LFPTOD(&first, res);
       return res;
}

l_fp
l_fp_init_from_double( double rhs)
{
       l_fp temp;
       DTOLFP(rhs, &temp);
       return temp;
}

void
l_fp_swap(l_fp * first, l_fp *second)
{
       l_fp temp = *second;

       *second = *first;
       *first = temp;

       return;
}

//----------------------------------------------------------------------
// testing the relational macros works better with proper predicate
// formatting functions; it slows down the tests a bit, but makes for
// readable failure messages.
//----------------------------------------------------------------------


bool
l_isgt (const l_fp first, const l_fp second)
{

       return L_ISGT(&first, &second);
}

bool
l_isgtu(const l_fp first, const l_fp second)
{

       return L_ISGTU(&first, &second);
}

bool
l_ishis(const l_fp first, const l_fp second)
{

       return L_ISHIS(&first, &second);
}

bool
l_isgeq(const l_fp first, const l_fp second)
{

       return L_ISGEQ(&first, &second);
}

bool
l_isequ(const l_fp first, const l_fp second)
{

       return L_ISEQU(&first, &second);
}


//----------------------------------------------------------------------
// test data table for add/sub and compare
//----------------------------------------------------------------------


static const lfp_hl addsub_tab[][3] = {
       // trivial idendity:
       {{0 ,0         }, { 0,0         }, { 0,0}},
       // with carry from fraction and sign change:
       {{-1,0x80000000}, { 0,0x80000000}, { 0,0}},
       // without carry from fraction
       {{ 1,0x40000000}, { 1,0x40000000}, { 2,0x80000000}},
       // with carry from fraction:
       {{ 1,0xC0000000}, { 1,0xC0000000}, { 3,0x80000000}},
       // with carry from fraction and sign change:
       {{0x7FFFFFFF, 0x7FFFFFFF}, {0x7FFFFFFF,0x7FFFFFFF}, {0xFFFFFFFE,0xFFFFFFFE}},
       // two tests w/o carry (used for l_fp<-->double):
       {{0x55555555,0xAAAAAAAA}, {0x11111111,0x11111111}, {0x66666666,0xBBBBBBBB}},
       {{0x55555555,0x55555555}, {0x11111111,0x11111111}, {0x66666666,0x66666666}},
       // wide-range test, triggers compare trouble
       {{0x80000000,0x00000001}, {0xFFFFFFFF,0xFFFFFFFE}, {0x7FFFFFFF,0xFFFFFFFF}}
};
static const size_t addsub_cnt = (sizeof(addsub_tab)/sizeof(addsub_tab[0]));
static const size_t addsub_tot = (sizeof(addsub_tab)/sizeof(addsub_tab[0][0]));



//----------------------------------------------------------------------
// epsilon estimation for the precision of a conversion double --> l_fp
//
// The error estimation limit is as follows:
//  * The 'l_fp' fixed point fraction has 32 bits precision, so we allow
//    for the LSB to toggle by clamping the epsilon to be at least 2^(-31)
//
//  * The double mantissa has a precsion 54 bits, so the other minimum is
//    dval * (2^(-53))
//
//  The maximum of those two boundaries is used for the check.
//
// Note: once there are more than 54 bits between the highest and lowest
// '1'-bit of the l_fp value, the roundtrip *will* create truncation
// errors. This is an inherent property caused by the 54-bit mantissa of
// the 'double' type.
double
eps(double d)
{

       return fmax(ldexp(1.0, -31), ldexp(fabs(d), -53));
}

//----------------------------------------------------------------------
// test addition
//----------------------------------------------------------------------
void
test_AdditionLR(void)
{
       size_t idx = 0;

       for (idx = 0; idx < addsub_cnt; ++idx) {
               l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
               l_fp op2 = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
               l_fp e_res = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
               l_fp res = l_fp_add(op1, op2);

               TEST_ASSERT_EQUAL_l_fp(e_res, res);
       }
       return;
}

void
test_AdditionRL(void)
{
       size_t idx = 0;

       for (idx = 0; idx < addsub_cnt; ++idx) {
               l_fp op2 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
               l_fp op1 = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
               l_fp e_res = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
               l_fp res = l_fp_add(op1, op2);

               TEST_ASSERT_EQUAL_l_fp(e_res, res);
       }
       return;
}


//----------------------------------------------------------------------
// test subtraction
//----------------------------------------------------------------------
void
test_SubtractionLR(void)
{
       size_t idx = 0;

       for (idx = 0; idx < addsub_cnt; ++idx) {
               l_fp op2 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
               l_fp e_res = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
               l_fp op1 = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
               l_fp res = l_fp_subtract(op1, op2);

               TEST_ASSERT_EQUAL_l_fp(e_res, res);
       }
       return;
}

void
test_SubtractionRL(void)
{
       size_t idx = 0;

       for (idx = 0; idx < addsub_cnt; ++idx) {
               l_fp e_res = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
               l_fp op2 = l_fp_init(addsub_tab[idx][1].h, addsub_tab[idx][1].l);
               l_fp op1 = l_fp_init(addsub_tab[idx][2].h, addsub_tab[idx][2].l);
               l_fp res = l_fp_subtract(op1, op2);

               TEST_ASSERT_EQUAL_l_fp(e_res, res);
       }
       return;
}

//----------------------------------------------------------------------
// test negation
//----------------------------------------------------------------------

void
test_Negation(void)
{
       size_t idx = 0;

       for (idx = 0; idx < addsub_cnt; ++idx) {
               l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
               l_fp op2 = l_fp_negate(op1);
               l_fp sum = l_fp_add(op1, op2);

               l_fp zero = l_fp_init(0, 0);

               TEST_ASSERT_EQUAL_l_fp(zero, sum);
       }
       return;
}



//----------------------------------------------------------------------
// test absolute value
//----------------------------------------------------------------------
void
test_Absolute(void)
{
       size_t idx = 0;

       for (idx = 0; idx < addsub_cnt; ++idx) {
               l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
               l_fp op2 = l_fp_abs(op1);

               TEST_ASSERT_TRUE(l_fp_signum(op2) >= 0);

               if (l_fp_signum(op1) >= 0)
                       op1 = l_fp_subtract(op1, op2);
               else
                       op1 = l_fp_add(op1, op2);

               l_fp zero = l_fp_init(0, 0);

               TEST_ASSERT_EQUAL_l_fp(zero, op1);
       }

       // There is one special case we have to check: the minimum
       // value cannot be negated, or, to be more precise, the
       // negation reproduces the original pattern.
       l_fp minVal = l_fp_init(0x80000000, 0x00000000);
       l_fp minAbs = l_fp_abs(minVal);
       TEST_ASSERT_EQUAL(-1, l_fp_signum(minVal));

       TEST_ASSERT_EQUAL_l_fp(minVal, minAbs);

       return;
}


//----------------------------------------------------------------------
// fp -> double -> fp rountrip test
//----------------------------------------------------------------------
void
test_FDF_RoundTrip(void)
{
       size_t idx = 0;

       // since a l_fp has 64 bits in it's mantissa and a double has
       // only 54 bits available (including the hidden '1') we have to
       // make a few concessions on the roundtrip precision. The 'eps()'
       // function makes an educated guess about the avilable precision
       // and checks the difference in the two 'l_fp' values against
       // that limit.

       for (idx = 0; idx < addsub_cnt; ++idx) {
               l_fp op1 = l_fp_init(addsub_tab[idx][0].h, addsub_tab[idx][0].l);
               double op2 = l_fp_convert_to_double(op1);
               l_fp op3 = l_fp_init_from_double(op2);

               l_fp temp = l_fp_subtract(op1, op3);
               double d = l_fp_convert_to_double(temp);
               TEST_ASSERT_DOUBLE_WITHIN(eps(op2), 0.0, fabs(d));
       }

       return;
}


//----------------------------------------------------------------------
// test the compare stuff
//
// This uses the local compare and checks if the operations using the
// macros in 'ntp_fp.h' produce mathing results.
// ----------------------------------------------------------------------
void
test_SignedRelOps(void)
{
       const lfp_hl * tv = (&addsub_tab[0][0]);
       size_t lc ;

       for (lc = addsub_tot - 1; lc; --lc, ++tv) {
               l_fp op1 = l_fp_init(tv[0].h, tv[0].l);
               l_fp op2 = l_fp_init(tv[1].h, tv[1].l);
               int cmp = l_fp_scmp(op1, op2);

               switch (cmp) {
               case -1:
                       //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
                       l_fp_swap(&op1, &op2);
                       //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
               case 1:
                       TEST_ASSERT_TRUE (l_isgt(op1, op2));
                       TEST_ASSERT_FALSE(l_isgt(op2, op1));

                       TEST_ASSERT_TRUE (l_isgeq(op1, op2));
                       TEST_ASSERT_FALSE(l_isgeq(op2, op1));

                       TEST_ASSERT_FALSE(l_isequ(op1, op2));
                       TEST_ASSERT_FALSE(l_isequ(op2, op1));
                       break;
               case 0:
                       TEST_ASSERT_FALSE(l_isgt(op1, op2));
                       TEST_ASSERT_FALSE(l_isgt(op2, op1));

                       TEST_ASSERT_TRUE (l_isgeq(op1, op2));
                       TEST_ASSERT_TRUE (l_isgeq(op2, op1));

                       TEST_ASSERT_TRUE (l_isequ(op1, op2));
                       TEST_ASSERT_TRUE (l_isequ(op2, op1));
                       break;
               default:
                       TEST_FAIL_MESSAGE("unexpected UCMP result: ");
               }
       }

       return;
}

void
test_UnsignedRelOps(void)
{
       const lfp_hl * tv =(&addsub_tab[0][0]);
       size_t lc;

       for (lc = addsub_tot - 1; lc; --lc, ++tv) {
               l_fp op1 = l_fp_init(tv[0].h, tv[0].l);
               l_fp op2 = l_fp_init(tv[1].h, tv[1].l);
               int cmp = l_fp_ucmp(op1, op2);

               switch (cmp) {
               case -1:
                       //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
                       l_fp_swap(&op1, &op2);
                       //printf("op1:%d %d, op2:%d %d\n",op1.l_uf,op1.l_ui,op2.l_uf,op2.l_ui);
               case 1:
                       TEST_ASSERT_TRUE (l_isgtu(op1, op2));
                       TEST_ASSERT_FALSE(l_isgtu(op2, op1));

                       TEST_ASSERT_TRUE (l_ishis(op1, op2));
                       TEST_ASSERT_FALSE(l_ishis(op2, op1));
                       break;
               case 0:
                       TEST_ASSERT_FALSE(l_isgtu(op1, op2));
                       TEST_ASSERT_FALSE(l_isgtu(op2, op1));

                       TEST_ASSERT_TRUE (l_ishis(op1, op2));
                       TEST_ASSERT_TRUE (l_ishis(op2, op1));
                       break;
               default:
                       TEST_FAIL_MESSAGE("unexpected UCMP result: ");
               }
       }

       return;
}

/*
*/

//----------------------------------------------------------------------
// that's all folks... but feel free to add things!
//----------------------------------------------------------------------

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.