Commits

ssnyder  committed b9d576f

Add FloatPacker.

  • Participants
  • Parent commits c41557b
  • Tags CxxUtils-00-00-43

Comments (0)

Files changed (6)

+2009-12-09  scott snyder  <snyder@bnl.gov>
+
+	* Tagging CxxUtils-00-00-42.
+	* CxxUtils/FloatPacker.h, src/FloatPacker.cxx,
+	test/FloatPacker_test.cxx, share/FloatPacker_test.ref: (new)
+
 2009-11-19  scott snyder  <snyder@bnl.gov>
 
 	* Tagging CxxUtils-00-00-41.

File CxxUtils/FloatPacker.h

+// This file's extension implies that it's C, but it's really -*- C++ -*-.
+// $Id$
+/**
+ * @file CxxUtils/FloatPacker.h
+ * @author scott snyder <snyder@bnl.gov>
+ * @date Nov, 2009, from earlier code.
+ * @brief Pack/unpack floating-point data from/to a given number of bits.
+ */
+
+
+#ifndef CXXUTILS_FLOATPACKER_H
+#define CXXUTILS_FLOATPACKER_H
+
+
+#include <string>
+#include <stdint.h>
+
+
+
+namespace CxxUtils {
+
+
+/**
+ * @brief Pack/unpack floating-point data from/to a given number of bits.
+ *
+ * The format is specified by the following parameters.
+ *
+ *   nbits - The total number of bits in the representation.
+ *   scale - Scale factor to apply before storing.
+ *   nmantissa - The number of bits to use for the mantissa and sign bit.
+ *   is_signed - Flag to tell if we should use a sign bit.
+ *   round - Flag to tell if we should round or truncate.
+ *
+ * From these we derive:
+ *
+ *   npack = nmantissa, if is_signed is false.
+ *         = nmantissa-1 if is_signed is true.
+ *   nexp = nbits - nmantissa
+ *
+ * The format consists of, in order from high bits to low bits:
+ *   - A sign bit, if is_signed is true.
+ *   - nexp bits of exponent information.
+ *   - npack bits of mantissa.
+ *
+ * The number is stored in normalized form, with an exponent bias
+ * of 2^(nexp-1).  But if the (biased) exponent is zero, then the
+ * mantissa is stored in denormalized form.  If nexp==0, this
+ * gives a fixed-point representation in the range [0,1).
+ * 0 is represented by all bits 0; if we have a sign bit, we can also
+ * represent -0 by all bits 0 except for the sign bit.
+ */
+class FloatPacker
+{
+public:
+  /// Type into which we pack.
+  typedef uint32_t Packdest;
+
+
+  /**
+   * @brief Constructor.
+   * @param nbits The number of bits in the packed representation.
+   * @param nmantissa The number of bits to use for the mantissa
+   *                  and sign bit.
+   * @param scale Divide the input number by this before packing.
+   * @param is_signed If true, then one mantissa bit is used for a sign.
+   * @param round If true, numbers will be rounded.
+   *              Otherwise, they will be truncated.
+   */
+  FloatPacker (int nbits,
+               int nmantissa,
+               double scale = 1,
+               bool is_signed = true,
+               bool round = false);
+
+
+  /**
+   * @brief Check to see if an error occurred.
+   * @param err[out] If an error occured, a description of it.
+   * @return True if an error occurred since the last call to @c errcheck.
+   */
+  bool errcheck (std::string& err) const;
+
+
+  /**
+   * @brief Pack a value.
+   * @param src Value to pack.
+   * @return The packed value.
+   *
+   * For now, we convert floats to doubles before packing.
+   */
+  Packdest pack (double src) const;
+
+
+  /**
+   * @brief Unpack the value @c VAL.
+   * @param val The packed data.  It should start with the low bit,
+   *            and any extraneous bits should have been masked off.
+   */
+  double unpack (Packdest val) const;
+
+
+private:
+  /// Number of bits in the mantissa + sign bit.
+  int m_nmantissa;
+
+  /// Scale factor for stored numbers.
+  double m_scale;
+
+  /// Should we use a sign bit?
+  bool m_is_signed;
+
+  /// Should we round instead of truncating?
+  bool m_round;
+
+  /// Number of bits in mantissa (exclusive of any sign bit).
+  int m_npack;
+
+  /// Mask with that many low bits set.
+  Packdest m_npack_ones;
+
+  /// Mask containing the sign bit (or 0 if there's no sign bit).
+  Packdest m_signmask;
+
+  /// Number of exponent bits.
+  int m_nexp;
+
+  /// Mask with that many low bits set.
+  Packdest m_nexp_ones;
+
+  /// Minimum exponent value.
+  int m_min_exp;
+
+  /// Maximum exponent value.
+  int m_max_exp;
+
+  /// Description of the last error.
+  mutable std::string m_lasterr;
+};
+
+
+} // namespace CxxUtils
+
+
+#endif // not CXXUTILS_FLOATPACKER_H

File cmt/requirements

 apply_pattern UnitTest_run unit_test=exctrace1
 apply_pattern UnitTest_run unit_test=exctrace2
 apply_pattern UnitTest_run unit_test=pointer_list
+apply_pattern UnitTest_run unit_test=FloatPacker
 
 end_private

File share/FloatPacker_test.ref

Empty file added.

File src/FloatPacker.cxx

+// $Id$
+/**
+ * @file CxxUtils/src/FloatPacker.cxx
+ * @author scott snyder <snyder@bnl.gov>
+ * @date Nov, 2009, from earlier code.
+ * @brief Pack/unpack floating-point data from/to a given number of bits.
+ */
+
+#include "CxxUtils/FloatPacker.h"
+#include <limits>
+#include <string>
+#include <sstream>
+#include <iomanip>
+#include <ieee754.h>  // ??? Is this standardized?
+#include <stdint.h>
+
+
+namespace {
+
+
+
+/**
+ * @brief We sometimes need to be able to access a double
+ *        as ints.
+ */
+union double_or_int {
+  ieee754_double d;
+  uint32_t i[2];
+};
+
+
+const int ieee754_double_bias = 0x3ff;
+const int ieee754_double_exponent_bits = 11;
+const int ieee754_double_mantissa0_bits = 20;
+const int ieee754_double_mantissa1_bits = 32;
+
+
+/// Abbreviation.
+typedef CxxUtils::FloatPacker::Packdest Packdest;
+
+/// Size of @c Packdest in bits.
+const int packdest_bits = std::numeric_limits<Packdest>::digits;
+
+
+// Handy constant: A Packdest with 1 in the high bit.
+const Packdest high_packdest_bit = (1U << (packdest_bits - 1));
+
+const Packdest ieee754_double_exponent_mask =
+  (1U << ieee754_double_exponent_bits) - 1;
+
+
+/**
+ * @brief  Return the largest (signed) integer that can be represented
+ *         With @c nbits bits.
+ *         Boundary case: nbits=0 should return -1.
+ * @param nbits Number of bits.
+ */
+inline
+int max_int (int nbits)
+{
+  return ((1U << nbits) >> 1) - 1;
+}
+
+
+/**
+ * @brief  Return the smallest (signed) integer that can be represented
+ *         With @c nbits bits.
+ *         Boundary case: nbits=0 should return -1.
+ * @param nbits Number of bits.
+ */
+inline
+int min_int (int nbits)
+{
+  return ((-1 << nbits) >> 1);
+}
+
+
+/**
+ * @brief Return a word with the low NBITS bits set.
+ * @param nbits Number of bits.
+ */
+Packdest ones (int nbits)
+{
+  if (nbits == packdest_bits)
+    return ~ static_cast<Packdest> (0);
+  return (static_cast<Packdest> (1) << nbits) - 1;
+}
+
+
+/**
+ * @brief Renormalize a denormal number.
+ * @param exponent[inout] The exponent of the number.
+ *                        Should initially be the denormal flag value.
+ *                        Will be modified in place.
+ * @param mantissa[inout] The mantissa of the number.
+ *                        Will be modified in place.
+ */
+void renormalize_denormal (int& exponent,
+                           Packdest& mantissa)
+{
+  if (mantissa == 0)
+    exponent -= packdest_bits; // Lost precision here.
+  else {
+    while ((mantissa & high_packdest_bit) == 0) {
+      --exponent;
+      mantissa <<= 1;
+    }
+    mantissa <<= 1;
+  }
+}
+
+
+/**
+ * @brief If the number given by @a exponent and @a mantissa is too small
+ *           to be represented by a normalized number in a representation
+ *           where @a min_exp is the exponent value flagging denormals, convert
+ *           it to a denormal representation.
+ *
+ * @param min_exp The denormal marker exponent value.
+ * @param round_bits Number of bits to consider for rounding.
+ * @param exponent[inout] The exponent of the number.
+ *                        Will be modified in place.
+ * param mantissa[inout] The mantissa of the number.
+ *                       Will be modified in place.
+ */
+void underflow_to_denormal (int min_exp,
+                            int round_bits,
+                            int& exponent,
+                            Packdest& mantissa)
+{
+  if (exponent <= min_exp) {
+    Packdest mantissa_in = mantissa;
+
+    // Denormalize the mantissa.
+    mantissa = (mantissa >> 1) | high_packdest_bit;
+
+    // Now shift it right.
+    int shift = min_exp - exponent;
+    if (shift < packdest_bits)
+      mantissa >>= shift;
+    else
+      mantissa = 0; // underflow to 0.
+
+    // Flag it as denormal.
+    exponent = min_exp;
+
+    // Handle rounding, if desired.
+    if (round_bits) {
+      //
+      //                                   +- packdest_bits - round_bits - 1
+      //                  |- round_bits  -|v
+      //    mantissa:     .................X....
+      //                       v- shift+1 -^
+      //    mantissa_in:  .....X................
+      //                       ^
+      //                       +- packdest_bits - round_bits + shift
+      //
+      int orig_pos = packdest_bits - round_bits + shift;
+      if (orig_pos < packdest_bits &&
+          ((static_cast<Packdest> (1) << orig_pos) & mantissa_in) != 0)
+      {
+        Packdest lsb = (static_cast<Packdest> (1) <<
+                        (packdest_bits - round_bits));
+        Packdest lsbmask = ~ (lsb - 1);
+
+        // ??? If we overflow here, it means we have to go back
+        //     to a normalized representation.  Just punt for now.
+        if ((mantissa & lsbmask) != lsbmask)
+          mantissa += lsb;
+      }
+    }
+  }
+}
+
+
+} // unnamed namespace
+
+
+namespace CxxUtils {
+
+
+/**
+ * @brief Constructor.
+ * @param nbits The number of bits in the packed representation.
+ * @param nmantissa The number of bits to use for the mantissa
+ *                  and sign bit.
+ * @param scale Divide the input number by this before packing.
+ * @param is_signed If true, then one mantissa bit is used for a sign.
+ * @param round If true, numbers will be rounded.
+ *              Otherwise, they will be truncated.
+ */
+FloatPacker::FloatPacker (int nbits,
+                          int nmantissa,
+                          double scale /*= 1*/,
+                          bool is_signed /*= true*/,
+                          bool round /*= false*/)
+  : m_nmantissa (nmantissa),
+    m_scale (scale),
+    m_is_signed (is_signed),
+    m_round (round)
+{
+  // scale==0 means not to scale.
+  // Use that instead of 1 since it's faster to test for 0.
+  if (scale == 1)
+    scale = 0;
+
+  // Set up other cached values.
+  m_npack = m_nmantissa;
+  if (m_is_signed)
+    --m_npack;
+
+  m_npack_ones = ones (m_npack);
+
+  // Sign bit mask.
+  if (m_is_signed)
+    m_signmask = 1U << (nbits - 1);
+  else
+    m_signmask = 0;
+
+  // Number of exponent bits.
+  m_nexp = nbits - m_nmantissa;
+  m_nexp_ones = ones (m_nexp);
+
+  // Minimum exponent value.
+  m_min_exp = min_int (m_nexp);
+
+  // Maximum exponent value.
+  m_max_exp = max_int (m_nexp);
+
+  if (m_npack < 1 || m_npack > nbits)
+    m_lasterr = "Bad number of mantissa bits.";
+}
+
+
+/**
+ * @brief Check to see if an error occurred.
+ * @param err[out] If an error occured, a description of it.
+ * @return True if an error occurred since the last call to @c errcheck.
+ */
+bool FloatPacker::errcheck (std::string& err) const
+{
+  if (!m_lasterr.empty()) {
+    err.swap (m_lasterr);
+    m_lasterr.clear();
+    return true;
+  }
+  return false;
+}
+
+
+/**
+ * @brief Pack a value.
+ * @param src Value to pack.
+ * @return The packed value.
+ *
+ * For now, we convert floats to doubles before packing.
+ */
+FloatPacker::Packdest FloatPacker::pack (double src) const
+{
+  double_or_int d;
+  d.d.d = src;
+
+  // Fast-path for zero.  (Purely an optimization.)
+  // Note: can't use a double compare here.  On some architectures (eg, MIPS)
+  // a denormal will compare equal to zero.
+  if (d.i[0] == 0 && d.i[1] == 0)
+    return 0;
+
+  // Check for NaN and infinity.
+  if (d.d.ieee.exponent == ieee754_double_exponent_mask) {
+    std::ostringstream os;
+    os << "Bad float number: " << src << " (" << std::setbase(16) << d.i[0]
+       << " " << d.i[1] << ")";
+    m_lasterr = os.str();
+    d.d.d = 0;
+  }
+
+  if (m_scale)
+    d.d.d /= m_scale;
+
+  bool was_negative = false;
+  if (d.d.ieee.negative != 0) {
+    if (m_is_signed) {
+      was_negative = true;
+      d.d.d = -d.d.d;
+    }
+    else {
+      // Don't complain on -0.
+      if (d.d.d < 0) {
+        std::ostringstream os;
+        os << "Float overflow during packing: " << src;
+        m_lasterr = os.str();
+      }
+      d.d.d = 0;
+    }
+  }
+
+  // Check for zero again.
+  // (Also need to preserve the sign; the scale division may
+  // have underflowed.)
+  if (d.i[0] == 0 && d.i[1] == 0) {
+    if (was_negative)
+      return m_signmask;
+    else
+      return 0;
+  }
+
+  // Get packdest_bits bits of mantissa.
+
+  Packdest mantissa =
+    (d.d.ieee.mantissa0 << (packdest_bits -
+                          ieee754_double_mantissa0_bits)) |
+    (d.d.ieee.mantissa1 >>
+     (ieee754_double_mantissa1_bits -
+      (packdest_bits - ieee754_double_mantissa0_bits)));
+
+  // Get the unbiased exponent.
+  int exponent =
+    static_cast<int> (d.d.ieee.exponent) - ieee754_double_bias;
+
+  // Do rounding, if requested.
+  if (m_round) {
+    Packdest lsbmask = (1 << (packdest_bits - m_npack));
+    int roundbit;
+    Packdest roundmask;
+    if (lsbmask > 1) {
+      roundbit = (mantissa & (lsbmask >> 1));
+      roundmask = ~ static_cast<Packdest> (roundbit - 1);
+    }
+    else {
+      roundbit = (d.d.ieee.mantissa1 &
+                  ((1 << ((ieee754_double_mantissa1_bits -
+                           (packdest_bits -
+                            ieee754_double_mantissa0_bits)) - 1))));
+      roundmask = ~ static_cast<Packdest> (0);
+    }
+
+    if (roundbit != 0) {
+      // Handle the case where it would overflow.
+      if ((mantissa & roundmask) == roundmask) {
+        mantissa >>= 1;
+        mantissa |= roundmask;
+        exponent += 1;
+      }
+
+      mantissa += lsbmask;
+    }
+  }
+
+  // If the number is too large, bitch, and reset to the largest number.
+  if (exponent > m_max_exp) {
+    std::ostringstream os;
+    os << "Float overflow during packing: " << src;
+    m_lasterr = os.str();
+    exponent = m_max_exp;
+    mantissa = static_cast<Packdest> (~0);
+  }
+
+  // Handle denormals.  (We've already handled the zero case.)
+  if (exponent == - ieee754_double_bias)
+    renormalize_denormal (exponent, mantissa);
+
+  // If the number is too small, denormalize, or underflow to 0.
+  underflow_to_denormal (m_min_exp, m_round ? m_npack: 0, exponent, mantissa);
+
+  // Pack in the mantissa bits.
+  Packdest dest = mantissa >> (packdest_bits - m_npack);
+
+  // The exponent, if desired.
+  if (m_nexp > 0)
+    dest |= ((exponent - m_min_exp) << m_npack);
+
+  // And the optional sign bit.
+  if (was_negative)
+    dest |= m_signmask;
+
+  return dest;
+}
+
+
+/**
+ * @brief Unpack the value @c VAL.
+ * @param val The packed data.  It should start with the low bit,
+ *            and any extraneous bits should have been masked off.
+ */
+double FloatPacker::unpack (Packdest val) const
+{
+  // Fast-path for 0.
+  if (val == 0)
+    return 0;
+
+  // Break apart the packed value.
+  bool was_negative = false;
+  if ((val & m_signmask) != 0)
+    was_negative = true;
+
+  double d;
+
+  // Fast path for fixed-point representations.
+  if (m_nexp == 0) {
+    Packdest mantissa = (val & m_npack_ones);
+    d = mantissa / ((double)m_npack_ones + 1);
+    if (was_negative)
+      d *= -1;
+  }
+  else {
+    // Get the mantissa.
+    Packdest mantissa = (val & m_npack_ones) << (packdest_bits - m_npack);
+
+    // General case.
+    // Get the exponent.
+    int exponent = ((val >> m_npack) & m_nexp_ones);
+    exponent += m_min_exp; // unbias.
+
+    ieee754_double dd;
+
+    // Handle denormals.
+    if (exponent == m_min_exp) {
+      // Maybe it was -0?
+      if (mantissa == 0) {
+        dd.d = 0;
+        if (was_negative)
+          dd.ieee.negative = 1;
+        return dd.d;
+      }
+
+      renormalize_denormal (exponent, mantissa);
+    }
+
+    // Complain about overflow.
+    if (exponent >= max_int (ieee754_double_exponent_bits)) {
+      std::ostringstream os;
+      os << "Overflow while unpacking float; exponent: " << exponent;
+      m_lasterr = os.str();
+      exponent = max_int (ieee754_double_exponent_bits) + 1;
+      mantissa = 0; // Infinity.
+    }
+
+    // Underflow into denormal.
+    underflow_to_denormal ( - ieee754_double_bias, 0,
+                            exponent, mantissa);
+
+    // Pack into a double.
+    dd.ieee.negative = was_negative ? 1 : 0;
+    dd.ieee.exponent = exponent + ieee754_double_bias;
+    dd.ieee.mantissa0 =
+      (mantissa >> (packdest_bits - ieee754_double_mantissa0_bits));
+    dd.ieee.mantissa1 =
+      (mantissa << (ieee754_double_mantissa0_bits -
+                    (packdest_bits - ieee754_double_mantissa1_bits)));
+    d = dd.d;
+  }
+
+  // Set the result.
+  if (m_scale)
+    d *= m_scale;
+  return d;
+}
+
+
+} // namespace CxxUtils

File test/FloatPacker_test.cxx

+// $Id$
+/**
+ * @file CxxUtils/test/FloatPacker_test.cxx
+ * @author scott snyder <snyder@bnl.gov>
+ * @date Nov, 2009, from earlier code.
+ * @brief Regression tests for FloatPacker.
+ */
+
+#include "CxxUtils/FloatPacker.h"
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include <ieee754.h>
+#include <string.h>
+
+
+using CxxUtils::FloatPacker;
+using std::cout;
+using std::abs;
+
+
+const int ieee754_double_exponent_bits = 11;
+
+
+bool bitwise_equal (double a, double b)
+{
+  return (memcmp (&a, &b, sizeof (a)) == 0);
+}
+
+
+bool almost_equal (double a, double b, int bits)
+{
+  return (abs ((a-b)/(a+b))) * (1<<bits) < 1;
+}
+
+
+void test1 ()
+{
+  // Set up some special numbers for testing.
+  ieee754_double d;
+  d.d = 0;
+  d.ieee.negative = 1;
+  double neg_zero = d.d;
+
+  d.d = 0;
+  d.ieee.exponent = (1<<ieee754_double_exponent_bits) - 1;
+  double infinity = d.d;
+
+#ifndef __alpha
+  const double denormal = 3e-320;
+  const double denormal2 = 3072e-320;
+#endif
+
+  FloatPacker tf1 (12, 12);
+
+  double unpacked;
+  FloatPacker::Packdest packed;
+  double out;
+  std::string err;
+
+  unpacked = 0.125;
+  packed = tf1.pack (unpacked);
+  assert (packed == 0x100);
+  out = tf1.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked = 0;
+  packed = tf1.pack (unpacked);
+  assert (packed == 0);
+  out = tf1.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked = -0.125;
+  packed = tf1.pack (unpacked);
+  assert (packed == 0x900);
+  out = tf1.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked = neg_zero;
+  packed = tf1.pack (unpacked);
+  assert (packed == 0x800);
+  out = tf1.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+  assert (!tf1.errcheck (err));
+
+  FloatPacker tf2 (28, 24, 2, false);
+
+  unpacked = 7.5;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0x9e00000);
+  out = tf2.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked = 7.5 + (0.5)/2/16/16/16/16/16;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0x9e00001);
+  out = tf2.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked = neg_zero;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0);
+  out = tf2.unpack (packed);
+  assert (out == 0);
+
+  unpacked = 0;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0);
+  out = tf2.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+  assert (!tf2.errcheck (err));
+
+  // This gets a warning.
+  unpacked = -1;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0);
+  out = tf2.unpack (packed);
+  assert (bitwise_equal (0, out));
+  assert (tf2.errcheck (err));
+  assert (err == "Float overflow during packing: -1");
+  assert (!tf2.errcheck (err));
+
+  // This gets a warning.
+  unpacked = infinity;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0);
+  out = tf2.unpack (packed);
+  assert (bitwise_equal (0, out));
+  assert (tf2.errcheck (err));
+  
+  assert (err == "Bad float number: inf (0 7ff00000)");
+
+  // This gets a warning.
+  unpacked = 100000;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0xfffffff);
+  out = tf2.unpack (packed);
+  assert (almost_equal (out, 512, 24));
+  assert (tf2.errcheck (err));
+  assert (err == "Float overflow during packing: 100000");
+
+  // Testing underflow to denormal.
+  unpacked = 0.5 / 32;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0x1000000);
+  out = tf2.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked /= 2;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0x0800000);
+  out = tf2.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked /= 2;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0x0400000);
+  out = tf2.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  // Underflow to zero.
+  unpacked = 1e-50;
+  packed = tf2.pack (unpacked);
+  assert (packed == 0);
+  out = tf2.unpack (packed);
+  assert (bitwise_equal (0, out));
+  assert (!tf2.errcheck (err));
+
+  FloatPacker tf3 (28, 12, 0, false);
+
+#ifndef __alpha
+  // Having problems with denormals on alpha.  Just skip for now.
+  unpacked = denormal;
+  packed = tf3.pack (unpacked);
+  assert (packed == 0x7be1000);
+  out = tf3.unpack (packed);
+  assert (bitwise_equal (0, out));
+
+  unpacked = denormal2;
+  packed = tf3.pack (unpacked);
+  assert (packed == 0x7be4400);
+  out = tf3.unpack (packed);
+#ifdef MIPS
+  // on mips, operating on denormals bashes them to zero.
+  assert (out == 0);
+#else
+  assert (almost_equal (out, unpacked, 3));
+#endif
+#endif
+  assert (!tf3.errcheck (err));
+
+  // Testing rounding.
+  FloatPacker tf6 (12, 8, 4, false, true);
+
+  unpacked = 990./1024; // This doesn't get rounded.
+  packed = tf6.pack (unpacked);
+  assert (packed == 0x5ef);
+  out = tf6.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked = 991./1024; // This does get rounded.
+  packed = tf6.pack (unpacked);
+  assert (packed == 0x5f0);
+  out = tf6.unpack (packed);
+  double tmp = 992./1024;
+  assert (bitwise_equal (tmp, out));
+
+  unpacked = 1023./1024; // Test mantissa overflowing after rounding.
+  packed = tf6.pack (unpacked);
+  assert (packed == 0x600);
+  out = tf6.unpack (packed);
+  tmp = 1.0;
+  assert (bitwise_equal (tmp, out));
+  assert (!tf6.errcheck (err));
+
+  // Test an unpacking problem with npack==packdest_bits.
+  FloatPacker tf7 (32, 32, 0, false);
+  unpacked = 1. / 65536 / 65536;
+  packed = tf7.pack (unpacked);
+  assert (packed == 0x1);
+  out = tf7.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+  assert (!tf7.errcheck (err));
+
+  // Test rounding during underflow_to_denormal.
+  FloatPacker tf8 (8, 8, 0, false, true);
+  unpacked = 0.5 + 1. / 256;
+  packed = tf8.pack (unpacked);
+  assert (packed == 0x81);
+  out = tf8.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked = 0.5 + 1. / 512;
+  packed = tf8.pack (unpacked);
+  assert (packed == 0x81);
+  out = tf8.unpack (packed);
+  tmp = 0.5 + 1. / 256;
+  assert (bitwise_equal (tmp, out));
+  assert (!tf8.errcheck (err));
+
+  // This gets a warning.
+  unpacked = 1023. / 1024;
+  packed = tf8.pack (unpacked);
+  assert (packed == 0xff);
+  d.d = unpacked;
+  assert (tf8.errcheck (err));
+  assert (err == "Float overflow during packing: 0.999023");  
+  out = tf8.unpack (packed);
+  tmp = 255. / 256;
+  assert (bitwise_equal (tmp, out));
+  assert (!tf8.errcheck (err));
+
+  // Test using bit 33 to control rounding.
+  FloatPacker tf9 (32, 32, 0, false, true);
+  unpacked = 0.5 + 1. / 65536 / 65536;
+  packed = tf9.pack (unpacked);
+  assert (packed == 0x80000001);
+  out = tf9.unpack (packed);
+  assert (bitwise_equal (unpacked, out));
+
+  unpacked = 0.5 + 1. / 65536 / 65536 / 2;
+  packed = tf9.pack (unpacked);
+  assert (packed == 0x80000001);
+  out = tf9.unpack (packed);
+  tmp = 0.5 + 1. / 65536 / 65536;
+  assert (bitwise_equal (tmp, out));
+  assert (!tf9.errcheck (err));
+
+  // This gets a warning.
+  packed = 0xf654000;
+  out = tf3.unpack (packed);
+  assert (tf3.errcheck (err));
+  assert (err == "Overflow while unpacking float; exponent: 30292");
+  assert (isinf (out));
+}
+
+
+int main ()
+{
+  test1 ();
+  return 0;
+}
+