X-Git-Url: https://git.ucc.asn.au/?p=ipdf%2Fcode.git;a=blobdiff_plain;f=src%2Ftests%2Frepresent.cpp;h=e5abef5b0e9524b265390b5d860132718ebcef35;hp=2fb524ad69d9a6be3ce9f7952225341ae216bc2d;hb=bf0f1cab25fa460ad4716ca60c11abab78a52cbd;hpb=805dd6babc18f4e16c297f0c20487418d5aa6bd8 diff --git a/src/tests/represent.cpp b/src/tests/represent.cpp index 2fb524a..e5abef5 100644 --- a/src/tests/represent.cpp +++ b/src/tests/represent.cpp @@ -3,6 +3,7 @@ #include #include #include +#include using namespace std; using namespace IPDF; @@ -11,6 +12,7 @@ using namespace IPDF; * This "tester" lets us construct Reals out of Floating Points with arbitrary sizes (but must be less than 64 bits currently) * Eventually it will also go the other way * This is originally for conceptual understanding but will ultimately become useful later on. + * Bahahaha */ @@ -26,73 +28,209 @@ long CountBits(long n) } /** - * Converts data represented as an (se, e, m, sm) floating point number to a Real. - * se = sign of exponent (1 bit) - * e = exponent (EMAX bits) - * m = mantissa (PREC bits) - * sm = sign of mantissa (1 bit) - * Total 1+EMAX+PREC+1 must not exceed 64 (at least, for the moment) + * Converts data represented as an (s, e, m) floating point number to a Real. + * Now actually correct for floats bigger than 8 bits! + * Now with 100% less memcpy and 100% more pointer arithmetic! + * Now with IEEE style encodings! (How many flops does it take to represent a float this way? TOO MANY!) + * Still probably breaks if you are not using an x86 based processor + * s = sign + * e = exponent (IEEE ``offset'' encoding) + * m = mantissa (IEEE encoding) */ -template Real BitsToReal(void * data) +template Real BitsToReal(void * data) { if (PREC + EMAX > 62) // we need 1 bit for the sign of the exponent, 1 bit for the sign Fatal("Can't do more than 62 bits (asked for %d + %d = %d)", PREC, EMAX, PREC+EMAX); - - // memcpy needs a whole number of bytes - // This means we need to do tricky bit shifting to get rid of parts of the exp/mantissa that overlap - size_t ebytes = ceil(double(EMAX+1.0) / 8.0); - size_t mbytes = ceil(double(PREC+1.0)/8.0); - size_t moffset = floor(double(EMAX+1.0) / 8.0); // offset to start of mantissa (nearest byte) - //Debug("ebytes %d, mbytes %d, moffset %d", ebytes, mbytes, moffset); - - // Get the exponent and its sign - uint64_t exponent = 0L; - bool exp_sign = false; - - memcpy(&exponent, data, ebytes); // Copy exponent - //Debug("exponent + garbage: %x", exponent); - exp_sign = (1L << (8*ebytes-1)) & exponent; - exponent &= ~(1L << (8*ebytes-1)); // mask out the sign - //Debug("exponent - masked sign: %x", exponent); - exponent = exponent >> (PREC+1); // shift out the extra bits (part of mantissa) - assert(CountBits(exponent) <= EMAX); //TODO: Remove once sure it actually works //TODO: Need more reliable sanity checks probably - - // Get the mantissa and its sign + + size_t mbytes = ceil(double(PREC)/8.0); + size_t mtrunc = floor(double(PREC)/8.0); uint64_t mantissa = 0L; - bool sign = false; - - memcpy(&mantissa, ((uint8_t*)(data) + moffset), mbytes); // copy data - //Debug("mantissa + garbage: %x", mantissa); - sign = mantissa & (1L); // get sign - mantissa = mantissa >> 1; // discard sign - mantissa = (mantissa << (8*sizeof(mantissa) - PREC)) >> (8*sizeof(mantissa) - PREC); - assert(CountBits(mantissa) <= PREC); - - /* - Debug("EMAX %u, PREC %u, BASE %u", EMAX, PREC, BASE); - Debug("EXP: %x", exponent); - Debug("MANTISSA: %x", mantissa); - Debug("EXP_SIGN: %x", (uint8_t)exp_sign); - Debug("MANTISSA_SIGN: %x", (uint8_t)sign); - */ + uint8_t * c = (uint8_t*)(data); + for (size_t i = 0; i < mbytes; ++i) + { + mantissa |= *(c++) << (i*8); + } + mantissa &= ~(~0L << PREC); - Real Q = (exp_sign) ? pow(Real(BASE), -Real(exponent)) : pow(Real(BASE), Real(exponent)); - Real x = Real(mantissa) * Q; + size_t ebytes = ceil(double(EMAX)/8.0); + uint64_t exponent = *(c-1) & (~0L << (PREC-8*mtrunc)); + for (size_t i = 0; i < ebytes; ++i) + { + exponent |= *(c++) << (i+1)*8; + } + exponent = exponent >> (PREC-8*mtrunc); + bool sign = exponent & (1LL << EMAX); + exponent &= ~(~0L << EMAX); + + // Floats are defined by + // x = (-1)^s * m * B^{e} + + // 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) + // So to conform with IEEE can't just treat the mantissa as an unsigned integer + + // Calculate the actual mantissa from the IEEE style encoding of the mantissa (sigh) + // This is obviously horribly inefficient but I guess that encoding not designed for ... whatever it is I'm doing here + Real M = (exponent != 0); + for (size_t i = 1; i < PREC; ++i) + { + if ((mantissa >> PREC-i) & 1) + M += Real(1) / powl(BASE, i); + } + + + //Debug("Mantissa: %x = %lf", mantissa, Float(M)); + //Debug("Exponent: %x = %u => %lf", exponent, exponent, Float(Real(exponent + 1) - Real(1LL << (EMAX-1)))); + //Debug("Sign: %x = %u", sign, sign); + + + Real x = M * powl(Real(BASE), Real(exponent + 1)-Real(1LL << (EMAX-1))); return (sign) ? -x : x; } +/** + * Modifies fraction to contain the fractional part of x + * Returns position of the first '1' bit + */ +size_t Fraction(Real x, uint64_t & fraction) +{ + x = abs(x) - (uint64_t)(abs(x)); + fraction = 0LL; + size_t offset = 0; + for (size_t i = 0; i < 8*sizeof(fraction); ++i) + { + x *= 2; + if (x >= 1.0) + { + offset = (offset == 0) ? i : offset; + fraction |= (1LL << (8*sizeof(fraction)-1-i)); + x -= 1.0; + } + if (x <= 0.0) + break; + } + return fraction; +} + +/** + * Copy mantissa and exponent into data + */ +template void BitsFromRepr(bool sign, uint64_t mantissa, uint64_t exponent, void * data) +{ + uint8_t * c = (uint8_t*)(data); + size_t i = 0; + for (i = 0; i < PREC; ++i) + { + *c |= ((mantissa >> i) & 1) << (i%8); + if ((i+1) % 8 == 0) + { + ++c; + } + } + for (; i < PREC+EMAX; ++i) + { + *c |= ((exponent >> (i-PREC)) & 1) << (i%8); + if ((i+1) % 8 == 0) + { + ++c; + } + } + *c |= (sign << (i%8)); +} + + + + /** * Performs the inverse of BitsToReal + * Hopefully + * It chooses the normalised representation */ -template void BitsFromReal(const Real & x, void * data) +template void BitsFromReal(Real x, void * data) { - bool sign; - uint64_t mantissa; - uint64_t exponent; if (PREC + EMAX > 62) Fatal("Can't do more than 62 bits (asked for %d + %d = %d)", PREC, EMAX, PREC+EMAX); + + uint64_t integer = (uint64_t)(abs(x)); + uint64_t fraction = 0LL; + const size_t max_bits = 8*sizeof(uint64_t); + bool sign = (x < 0); + x = abs(x) - integer; - //TODO: Implement + size_t fraction_leading_zeroes = max_bits; + for (size_t i = 0; i < max_bits; ++i) + { + x *= 2; + if (x >= 1.0) + { + if (fraction_leading_zeroes > i) + fraction_leading_zeroes = i; + x -= 1; + fraction |= (1LL << (max_bits-1-i)); + } + if (x <= 0.0) + break; + } + + uint64_t mantissa = 0LL; + int64_t biased_exponent = 0LL; + + // Combine integer and fraction into mantissa + if (integer == 0) + { + //Debug("0. %lx (%d)", fraction, fraction_leading_zeroes); + // Mantissa must have leading '1'; shift to get leading 1, then shift back to mantissa position + mantissa = (fraction << fraction_leading_zeroes) >> (max_bits - PREC); + biased_exponent -= fraction_leading_zeroes; + // If the unbiased exponent would be negative use a denormalised number instead + if (biased_exponent < -(1LL << (EMAX-1))+1) + { + mantissa = mantissa >> (-(1LL << (EMAX-1))+1 - biased_exponent); + biased_exponent = -(1LL << (EMAX-1))+1; + } + } + else + { + + size_t integer_leading_zeroes = 0; + for (integer_leading_zeroes = 0; integer_leading_zeroes < max_bits; ++integer_leading_zeroes) + { + if ((integer >> (max_bits-1-integer_leading_zeroes)) & 1LL) + break; + } + //Debug("%.16lx . %.16lx (%li, %li)", integer,fraction, integer_leading_zeroes, fraction_leading_zeroes); + + mantissa = (integer << (integer_leading_zeroes)); // shift integer part + mantissa |= (fraction >> (max_bits-integer_leading_zeroes)); // append the fraction after it + //Debug("%.16lx (%.16lx << %i)", mantissa, fraction, (fraction_leading_zeroes - max_bits + integer_leading_zeroes)); + mantissa = mantissa >> (max_bits - PREC); // shift back to the mantissa position + biased_exponent = (max_bits - integer_leading_zeroes); // calculate exponent + } + + // If the unbiased exponent would be non-zero, the IEEE has the leading 1 as implied; remove it and decrease exponent to compensate + if (biased_exponent != -(1LL << (EMAX-1))+1) + { + mantissa &= ~(1LL << max_bits - PREC); + mantissa = mantissa << 1; + biased_exponent -= 1; + } + + uint64_t exponent = (uint64_t)(biased_exponent + (1LL << (EMAX-1))-1); // offset the exponent + + // Debug + /* + Real M = (exponent != 0); + for (size_t i = 1; i < PREC; ++i) + { + if ((mantissa >> PREC-i) & 1) + M += Real(1) / powl(BASE, i); + } + + Debug("Exponent: %lx %lu -> %li", exponent, exponent, biased_exponent); + Debug("Mantissa: %lx -> %f", mantissa, Float(M)); + Debug("Both: %lx", mantissa | (exponent << (PREC))); + */ + // Now copy the bits in + BitsFromRepr(sign, mantissa, exponent, data); } @@ -103,23 +241,28 @@ int main(int argc, char ** argv) printf("# Convert custom floats to a Real\n"); printf("# a\thex(a)\tReal(a)\tdelta(last2)\n"); - typedef pair Pear; + typedef pair Pear; list space; - - for (uint8_t a = 0x00; a < 0xFF; ++a) + uint16_t a0 = 0x0000; + for (uint16_t a = a0; a < 0xFFFF; ++a) { - Real x = BitsToReal<2,4>(&a); - space.push_back(pair(a, x)); + Real x = BitsToReal<5,10>(&a); + uint16_t b = 0; BitsFromReal<5,10>(x, &b); + Real y = BitsToReal<5,10>(&b); + if (y != x) + { + Fatal("%x -> %lf -> %x -> %lf", a, Float(x), b, Float(y)); + } + space.push_back(Pear(a, x)); } space.sort([=](const Pear & a, const Pear & b){return a.second < b.second;}); Real prev; for (list::iterator i = space.begin(); i != space.end(); ++i) { - printf("%u\t%x\t%f", i->first, i->first, Float(i->second)); + printf("%u\t%.2x\t%.20lf", i->first, i->first, Float(i->second)); if (i != space.begin()) - { - printf("\t%f\n", Float(i->second - prev)); - } + printf("\t%f", Float(i->second - prev)); + printf("\n"); prev = i->second; } }