src/tests/represent.cpp

   1 #include "main.h"
   2 #include "real.h"
   3 #include <cmath>
   4 #include <cassert>
   5 #include <list>
   6 #include <bitset>
   7
   8 using namespace std;
   9 using namespace IPDF;
  10
  11 /**
  12  * This "tester" lets us construct Reals out of Floating Points with arbitrary sizes (but must be less than 64 bits currently)
  13  * Eventually it will also go the other way
  14  * This is originally for conceptual understanding but will ultimately become useful later on.
  15  * Bahahaha
  16  */
  17
  18
  19 /**
  20  * From http://stackoverflow.com/questions/109023/how-to-count-the-number-of-set-bits-in-a-32-bit-integer
  21  */
  22 long CountBits(long n)
  23 {
  24         unsigned int c; // c accumulates the total bits set in v
  25         for (c = 0; n; c++)
  26                 n &= n - 1; // clear the least significant bit set
  27         return c;
  28 }
  29
  30 /**
  31  * Converts data represented as an (s, e, m) floating point number to a Real.
  32  * Now actually correct for floats bigger than 8 bits!
  33  * Now with 100% less memcpy and 100% more pointer arithmetic!
  34  * Now with IEEE style encodings! (How many flops does it take to represent a float this way? TOO MANY!)
  35  * Still probably breaks if you are not using an x86 based processor
  36  * s = sign
  37  * e = exponent (IEEE ``offset'' encoding)
  38  * m = mantissa (IEEE encoding)
  39  */
  40 template <size_t EMAX, size_t PREC, size_t BASE = 2> Real BitsToReal(void * data)
  41 {
  42         if (PREC + EMAX > 62) // we need 1 bit for the sign of the exponent, 1 bit for the sign
  43                 Fatal("Can't do more than 62 bits (asked for %d + %d = %d)", PREC, EMAX, PREC+EMAX);
  44
  45         size_t mbytes = ceil(double(PREC)/8.0);
  46         size_t mtrunc = floor(double(PREC)/8.0);
  47         uint64_t mantissa = 0L;
  48         uint8_t * c = (uint8_t*)(data);
  49         for (size_t i = 0; i < mbytes; ++i)
  50         {
  51                 mantissa |= *(c++) << (i*8);
  52         }
  53         mantissa &= ~(~0L << PREC);
  54
  55         size_t ebytes = ceil(double(EMAX)/8.0);
  56         uint64_t exponent = *(c-1) & (~0L << (PREC-8*mtrunc));
  57         for (size_t i = 0; i < ebytes; ++i)
  58         {
  59                 exponent |= *(c++) << (i+1)*8;
  60         }
  61         exponent = exponent >> (PREC-8*mtrunc);
  62         bool sign = exponent & (1LL << EMAX);
  63         exponent &= ~(~0L << EMAX);
  64
  65         // Floats are defined by
  66         // x = (-1)^s * m * B^{e}
  67
  68         // In IEEE things get a little crazy though, the mantissa is stored reversed and the digits are summed indexed using negative powers of 2 (not positive)
  69         // So to conform with IEEE can't just treat the mantissa as an unsigned integer
  70
  71         // Calculate the actual mantissa from the IEEE style encoding of the mantissa (sigh)
  72         // This is obviously horribly inefficient but I guess that encoding not designed for ... whatever it is I'm doing here
  73         Real M = (exponent != 0);
  74         for (size_t i = 1; i < PREC; ++i)
  75         {
  76                 if ((mantissa >> PREC-i) & 1)
  77                         M += Real(1) / powl(BASE, i);
  78         }
  79
  80
  81         //Debug("Mantissa: %x = %lf", mantissa, Float(M));
  82         //Debug("Exponent: %x = %u => %lf", exponent, exponent, Float(Real(exponent + 1) - Real(1LL << (EMAX-1))));
  83         //Debug("Sign: %x = %u", sign, sign);
  84
  85
  86         Real x = M * powl(Real(BASE), Real(exponent + 1)-Real(1LL << (EMAX-1)));
  87         return (sign) ? -x : x;
  88 }
  89
  90 /**
  91  * Modifies fraction to contain the fractional part of x
  92  * Returns position of the first '1' bit
  93  */
  94 size_t Fraction(Real x, uint64_t & fraction)
  95 {
  96         x = abs(x) - (uint64_t)(abs(x));
  97         fraction = 0LL;
  98         size_t offset = 0;
  99         for (size_t i = 0; i < 8*sizeof(fraction); ++i)
 100         {
 101                 x *= 2;
 102                 if (x >= 1.0)
 103                 {
 104                         offset = (offset == 0) ? i : offset;
 105                         fraction |= (1LL << (8*sizeof(fraction)-1-i));
 106                         x -= 1.0;
 107                 }
 108                 if (x <= 0.0)
 109                         break;
 110         }
 111         return fraction;
 112 }
 113
 114 /**
 115  * Copy mantissa and exponent into data
 116  */
 117 template <size_t EMAX, size_t PREC, size_t BASE = 2> void BitsFromRepr(bool sign, uint64_t mantissa, uint64_t exponent, void * data)
 118 {
 119         uint8_t * c = (uint8_t*)(data);
 120         size_t i = 0;
 121         for (i = 0; i < PREC; ++i)
 122         {
 123                 *c |= ((mantissa >> i) & 1) << (i%8);
 124                 if ((i+1) % 8 == 0)
 125                 {
 126                         ++c;
 127                 }
 128         }
 129         for (; i < PREC+EMAX; ++i)
 130         {
 131                 *c |= ((exponent >> (i-PREC)) & 1) << (i%8);
 132                 if ((i+1) % 8 == 0)
 133                 {
 134                         ++c;
 135                 }
 136         }
 137         *c |= (sign << (i%8));
 138 }
 139
 140
 141
 142
 143 /**
 144  * Performs the inverse of BitsToReal
 145  * Hopefully
 146  * It chooses the normalised representation
 147  */
 148 template <size_t EMAX, size_t PREC, size_t BASE = 2> void BitsFromReal(Real x, void * data)
 149 {
 150         if (PREC + EMAX > 62)
 151                 Fatal("Can't do more than 62 bits (asked for %d + %d = %d)", PREC, EMAX, PREC+EMAX);
 152
 153         uint64_t integer = (uint64_t)(abs(x));
 154         uint64_t fraction = 0LL;
 155         const size_t max_bits = 8*sizeof(uint64_t);
 156         bool sign = (x < 0);
 157         x = abs(x) - integer;
 158
 159         size_t fraction_leading_zeroes = max_bits;
 160         for (size_t i = 0; i < max_bits; ++i)
 161         {
 162                 x *= 2;
 163                 if (x >= 1.0)
 164                 {
 165                         if (fraction_leading_zeroes > i)
 166                                 fraction_leading_zeroes = i;
 167                         x -= 1;
 168                         fraction |= (1LL << (max_bits-1-i));
 169                 }
 170                 if (x <= 0.0)
 171                         break;
 172         }
 173
 174         uint64_t mantissa = 0LL;
 175         int64_t biased_exponent = 0LL;
 176
 177         // Combine integer and fraction into mantissa
 178         if (integer == 0)
 179         {
 180                 //Debug("0. %lx (%d)", fraction, fraction_leading_zeroes);
 181                 // Mantissa must have leading '1'; shift to get leading 1, then shift back to mantissa position
 182                 mantissa = (fraction << fraction_leading_zeroes) >> (max_bits - PREC);
 183                 biased_exponent -= fraction_leading_zeroes;
 184                 // If the unbiased exponent would be negative use a denormalised number instead
 185                 if (biased_exponent < -(1LL << (EMAX-1))+1)
 186                 {
 187                         mantissa = mantissa >> (-(1LL << (EMAX-1))+1 - biased_exponent);
 188                         biased_exponent = -(1LL << (EMAX-1))+1;
 189                 }
 190         }
 191         else
 192         {
 193
 194                 size_t integer_leading_zeroes = 0;
 195                 for (integer_leading_zeroes = 0; integer_leading_zeroes < max_bits; ++integer_leading_zeroes)
 196                 {
 197                         if ((integer >> (max_bits-1-integer_leading_zeroes)) & 1LL)
 198                                 break;
 199                 }
 200                 //Debug("%.16lx . %.16lx (%li, %li)", integer,fraction,  integer_leading_zeroes, fraction_leading_zeroes);
 201
 202                 mantissa = (integer << (integer_leading_zeroes)); // shift integer part
 203                 mantissa |= (fraction >> (max_bits-integer_leading_zeroes)); // append the fraction after it
 204                 //Debug("%.16lx (%.16lx << %i)", mantissa,  fraction, (fraction_leading_zeroes - max_bits + integer_leading_zeroes));
 205                 mantissa = mantissa >> (max_bits - PREC); // shift back to the mantissa position
 206                 biased_exponent = (max_bits - integer_leading_zeroes); // calculate exponent
 207         }
 208
 209         // If the unbiased exponent would be non-zero, the IEEE has the leading 1 as implied; remove it and decrease exponent to compensate
 210         if (biased_exponent != -(1LL << (EMAX-1))+1)
 211         {
 212                 mantissa &= ~(1LL << max_bits - PREC);
 213                 mantissa = mantissa << 1;
 214                 biased_exponent -= 1;
 215         }
 216
 217         uint64_t exponent = (uint64_t)(biased_exponent + (1LL << (EMAX-1))-1); // offset the exponent
 218
 219         // Debug
 220         /*
 221         Real M = (exponent != 0);
 222         for (size_t i = 1; i < PREC; ++i)
 223         {
 224                 if ((mantissa >> PREC-i) & 1)
 225                         M += Real(1) / powl(BASE, i);
 226         }
 227
 228         Debug("Exponent: %lx %lu -> %li", exponent, exponent, biased_exponent);
 229         Debug("Mantissa: %lx -> %f", mantissa, Float(M));
 230         Debug("Both: %lx", mantissa | (exponent << (PREC)));
 231         */
 232         // Now copy the bits in
 233         BitsFromRepr<EMAX, PREC, BASE>(sign, mantissa, exponent, data);
 234
 235 }
 236
 237
 238
 239 int main(int argc, char ** argv)
 240 {
 241         printf("# Convert custom floats to a Real\n");
 242         printf("# a\thex(a)\tReal(a)\tdelta(last2)\n");
 243
 244         typedef pair<uint16_t, Real> Pear;
 245         list<Pear> space;
 246         uint16_t a0 = 0x0000;
 247         for (uint16_t a = a0; a < 0xFFFF; ++a)
 248         {
 249                 Real x = BitsToReal<5,10>(&a);
 250                 uint16_t b = 0; BitsFromReal<5,10>(x, &b);
 251                 Real y = BitsToReal<5,10>(&b);
 252                 if (y != x)
 253                 {
 254                         Fatal("%x -> %lf -> %x -> %lf", a, Float(x), b, Float(y));
 255                 }
 256                 space.push_back(Pear(a, x));
 257         }
 258         space.sort([=](const Pear & a, const Pear & b){return a.second < b.second;});
 259         Real prev;
 260         for (list<Pear>::iterator i = space.begin(); i != space.end(); ++i)
 261         {
 262                 printf("%u\t%.2x\t%.20lf", i->first, i->first, Float(i->second));
 263                 if (i != space.begin())
 264                         printf("\t%f", Float(i->second - prev));
 265                 printf("\n");
 266                 prev = i->second;
 267         }
 268 }