+#include <cmath>
+#include "binfloat.h"
+
+using namespace std;
+
+
+namespace Msp {
+namespace DataFile {
+
+BinFloat BinFloat::explode(UInt64 value, const Bits &bits)
+{
+ UInt64 mantissa_mask = (UInt64(1)<<bits.mantissa)-1;
+ int exponent_mask = (1<<bits.exponent)-1;
+
+ BinFloat bf;
+ // Extract biased exponent and sign
+ bf.exponent = (value>>bits.mantissa)&exponent_mask;
+ bf.sign = value>>(bits.mantissa+bits.exponent);
+ bf.infinity = (bf.exponent==exponent_mask);
+
+ if(bf.exponent==0 || bf.infinity)
+ // Zeroes and infinities have zero mantissa
+ bf.mantissa = 0;
+ else
+ {
+ // Extract mantissa, add the implied one and align it to high bits
+ bf.mantissa = (value&mantissa_mask) | (UInt64(1)<<bits.mantissa);
+ bf.mantissa <<= 63-bits.mantissa;
+ }
+
+ // Unbias the exponent
+ bf.exponent -= exponent_mask>>1;
+
+ return bf;
+}
+
+UInt64 BinFloat::compose(const Bits &bits)
+{
+ UInt64 mantissa_mask = (UInt64(1)<<bits.mantissa)-1;
+ int exponent_mask = (1<<bits.exponent)-1;
+
+ int biased_exponent = exponent+(exponent_mask>>1);
+ // Shift down and round the mantissa
+ UInt64 rounded_mantissa = ((mantissa>>(62-bits.mantissa))+1)>>1;
+ // If the integer part is greater than one, we need to use a higher exponent
+ if((rounded_mantissa>>bits.mantissa)>1)
+ ++biased_exponent;
+
+ if(biased_exponent>=exponent_mask || infinity)
+ // Overflow, return infinity
+ return UInt64(sign<<bits.exponent | exponent_mask)<<bits.mantissa;
+ else if(biased_exponent<=0 || !mantissa)
+ // Underflow, return zero
+ return 0;
+ else
+ {
+ UInt64 value = rounded_mantissa&mantissa_mask;
+ value |= UInt64(biased_exponent)<<bits.mantissa;
+ value |= UInt64(sign)<<(bits.mantissa+bits.exponent);
+ return value;
+ }
+}
+
+
+// exponent = log_2(bits)*3-7
+BinFloat::Bits::Bits(unsigned bits):
+ exponent(log(bits)*4.3281-7),
+ mantissa(bits-exponent-1)
+{ }
+
+} // namespace DataFile
+} // namespace Msp