@@@ fltfmt wip

[mLib] / utils / fltfmt.c

/* -*-c-*-
 *
 * Floating-point format conversions
 *
 * (c) 2024 Straylight/Edgeware
 */

/*----- Licensing notice --------------------------------------------------*
 *
 * This file is part of the mLib utilities library.
 *
 * mLib is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Library General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * mLib is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with mLib.  If not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
 * USA.
 */

/*----- Header files ------------------------------------------------------*/

#include "config.h"

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

#include "alloc.h"
#include "arena.h"
#include "bits.h"
#include "fltfmt.h"
#include "growbuf.h"
#include "maths.h"

/*----- Preliminary hacking -----------------------------------------------*/

/* The native-format conversions are -- at least if the format is
 * unrecognized -- dependent on the implementation's rounding.  Our own
 * rounding mode specifications don't fit into the framework very well, but I
 * still want to respect the prevailing rounding mode.
 *
 * The `proper' way to do this is with %|#pragma STDC FENV_ACCESS|%.  But
 * that doesn't actually work on GCC, or on Clang from not too long ago.  So
 * use compiler-specific hacking to support this.
 */
#if GCC_VERSION_P(4, 4)
#  pragma GCC optimize "-frounding-math"
#elif CLANG_VERSION_P(11, 0) && !CLANG_VERSION_P(12, 0)
#  pragma clang optimize "-frounding-math"
#elif GCC_VERSION_P(0, 0) || \
	(CLANG_VERSION_P(0, 0) && !CLANG_VERSION_P(12, 0))
   /* We just lose. */
#elif __STDC_VERSION__ >= 199001
#  include <fenv.h>
#  pragma STDC FENV_ACCESS ON
#endif

/*----- Some useful constants ---------------------------------------------*/

#define B31 0x80000000u			/* just bit 31 */
#define B30 0x40000000u			/* just bit 30 */

#define SH32 4294967296.0		/* 2^32 as floating-point */

/*----- Useful macros -----------------------------------------------------*/

#define B32(k) ((uint32)1 << (k))
#define M32(k) (B32(k) - 1)

#define FINITEP(x) (!((x)->f&(FLTF_INF | FLTF_NANMASK | FLTF_ZERO)))

/*----- Utility functions -------------------------------------------------*/

/* --- @clz32@ --- *
 *
 * Arguments:	@uint32 x@ = a 32-bit value
 *
 * Returns:	The number of leading zeros in @x@, i.e., the nonnegative
 *		integer %$n$% such that %$2^{31} \le 2^n x < 2^{32}$%.
 *		Returns a nonsensical value if %$x = 0$%.
 */

static unsigned clz32(uint32 x)
{
  unsigned n = 0;

  /* Divide and conquer.  If the top half of the bits are clear, then there
   * must be at least 16 leading zero bits, so accumulate and shift.  Repeat
   * for smaller powers of two.
   *
   * This ends up returning 31 if %$x = 0$%, but don't rely on this.
   */
  if (!(x&0xffff0000)) { x <<= 16; n += 16; }
  if (!(x&0xff000000)) { x <<=  8; n +=  8; }
  if (!(x&0xf0000000)) { x <<=  4; n +=  4; }
  if (!(x&0xc0000000)) { x <<=  2; n +=  2; }
  if (!(x&0x80000000)) {	   n +=  1; }
  return (n);
}

/* --- @ctz32@ --- *
 *
 * Arguments:	@uint32 x@ = a 32-bit value
 *
 * Returns:	The number of trailing zeros in @x@, i.e., the nonnegative
 *		integer %$n$% such that %$x/2^n$% is an odd integer.
 *		Returns a nonsensical value if %$x = 0$%.
 */

static unsigned ctz32(uint32 x)
{
#ifdef CTZ_TRADITIONAL

  unsigned n = 0;

  /* Divide and conquer.  If the bottom half of the bits are clear, then
   * there must be at least 16 trailing zero bits, so accumulate and shift.
   * Repeat for smaller powers of two.
   *
   * This ends up returning 31 if %$x = 0$%, but don't rely on this.
   */
  if (!(x&0x0000ffff)) { x >>= 16; n += 16; }
  if (!(x&0x000000ff)) { x >>=  8; n +=  8; }
  if (!(x&0x0000000f)) { x >>=  4; n +=  4; }
  if (!(x&0x00000003)) { x >>=  2; n +=  2; }
  if (!(x&0x00000001)) {	   n +=  1; }
  return (n);

#else

  static unsigned char tab[] =
    /*
       ;;; Compute the decoding table for the de Bruijn sequence trick below.

       (let ((db #x04653adf)
	     (rv (make-vector 32 nil)))
	 (dotimes (i 32)
	   (aset rv (logand (ash db (- i 27)) #x1f) i))
	 (save-excursion
	   (goto-char (point-min))
	   (search-forward (concat "***" "BEGIN ctz32tab" "***"))
	   (beginning-of-line 2)
	   (delete-region (point)
			  (progn
			    (search-forward "***END***")
			    (beginning-of-line)
			    (point)))
	   (dotimes (i 32)
	     (cond ((zerop i) (insert "    { "))
		   ((zerop (mod i 16)) (insert ",\n      "))
		   ((zerop (mod i 4)) (insert ",  "))
		   (t (insert ", ")))
	     (insert (format "%2d" (aref rv i))))
	   (insert " };\n")))
     */

    /* ***BEGIN ctz32tab*** */
    {  0,  1,  2,  6,   3, 11,  7, 16,   4, 14, 12, 21,   8, 23, 17, 26,
      31,  5, 10, 15,  13, 20, 22, 25,  30,  9, 19, 24,  29, 18, 28, 27 };
    /* ***END*** */

  /* Sneaky trick.  Two's complement negation (which you effectively get
   * using C unsigned arithmetic, whether you like it or not) complements all
   * of the bits of an operand more significant than the least significant
   * set bit.  Therefore, this bit is the only one set in both %$x$% and
   * %$-x$%, so @x&-x@ will isolate it for us.
   *
   * The magic number @0x04653adf@ is a %%\emph{de Bruijn} number%%: every
   * group of five consecutive bits is distinct, including groups which `wrap
   * around', including some low bits and some high bits.  Multiplying this
   * number by a power of two is equivalent to a left shift; and, because the
   * top five bits are all zero, the most significant five bits of the
   * product are the same as if we'd applied a rotation.  The result is that
   * we end up with a distinctive pattern in those bits which perfectly
   * diagnose each shift from 0 up to 31, which we can decode using a table.
   *
   * David Seal described a similar trick -- using the six-bit pattern
   * generated by the constant @0x0450fbaf@ -- in `comp.sys.arm' in 1994;
   * this constant was particularly convenient to multiply by on early ARM
   * processors.  The use of a de Bruijn number is described in Henry
   * Warren's %%\emph{Hacker's Delight}%%.
   */
  return (tab[((x&-x)*0x04653adf >> 27)&0x1f]);

#endif
}

/* --- @shl@, @shr@ --- *
 *
 * Arguments:	@uint32 *z@ = destination vector
 *		@const uint32 *x@ = source vector
 *		@size_t sz@ = size of source vector, in elements
 *		@unsigned n@ = number of bits to shift by; must be less than
 *			32
 *
 * Returns:	The bits shifted off the end of the vector.
 *
 * Use:		Shift a vector of 32-bit words left (@shl@) or right (@shr@)
 *		by some number of bit positions.  These functions work
 *		correctly if @z@ and @x@ are the same pointer, but not if
 *		they otherwise overlap.
 */

static uint32 shl(uint32 *z, const uint32 *x, size_t sz, unsigned n)
{
  size_t i;
  uint32 t, u;
  unsigned r;

  if (!n) {
    for (i = 0; i < sz; i++) z[i] = x[i];
    t = 0;
  } else {
    r = 32 - n;
    for (t = 0, i = sz; i--; )
      { u = x[i]; z[i] = ((u << n) | t)&MASK32; t = u >> r; }
  }
  return (t);
}

static uint32 shr(uint32 *z, const uint32 *x, size_t sz, unsigned n)
{
  size_t i;
  uint32 t, u;
  unsigned r;

  if (!n) {
    for (i = 0; i < sz; i++) z[i] = x[i];
    t = 0;
  } else {
    r = 32 - n;
    for (t = 0, i = 0; i < sz; i++)
      { u = x[i]; z[i] = ((u >> n) | t)&MASK32; t = u << r; }
  }
  return (t);
}

/* --- @sigbits@ --- *
 *
 * Arguments:	@const struct floatbits *x@ = decoded floating-point number
 *
 * Returns:	The number of significant digits in @x@'s fraction.  This
 *		will be zero if @x@ is zero or infinite.
 */

static unsigned sigbits(const struct floatbits *x)
{
  unsigned i;
  uint32 w;

  if (x->f&(FLTF_ZERO | FLTF_INF)) return (0);
  i = x->n;
  for (;;) {
    if (!i) return (0);
    w = x->frac[--i]; if (w) return (32*(i + 1) - ctz32(w));
  }
}

/* --- @ms_set_bit@ --- *
 *
 * Arguments:	@const uint32 *x@ = pointer to the %%\emph{end}%% of the
 *			buffer
 *		@unsigned from, to@ = lower (inclusive) and upper (exclusive)
 *			bounds on the region of bits to inspect
 *
 * Returns:	Index of the most significant set bit, or @ALLCLEAR@.
 *
 * Use:		For the (rather unusual) purposes of this function, the bits
 *		of the input are numbered from zero, being the least
 *		significant bit of @x[-1]@, upwards through more significant
 *		bits of @x[-1]@, and then through @x[-2]@ and so on.
 *
 *		If all of the bits in the half-open region are clear then
 *		@ALLCLEAR@ is returned; otherwise, the return value is the
 *		index of the most significant set bit in the region.  Note
 *		that @ALLCLEAR@ is equal to @UINT_MAX@: since this is the
 *		largest possible value of @to@, and the upper bound is
 *		exclusive, this cannot be the index of a bit in the region.
 */

#define ALLCLEAR UINT_MAX
static unsigned ms_set_bit(const uint32 *x, unsigned from, unsigned to)
{
  unsigned n0, n, b, base;
  uint32 m, w;

  /*		    <--- increasing indices <---
   *
   *     ---+-------+-------+-------+-------+-------+---
   * ...S   |///|   |       |       |       |    |//|   S...
   *     ---+-------+-------+-------+-------+-------+---
   */

  /* If the region is empty then it's technically true that all of the bits
   * are zero.  It's important to be able to do answer the case where
   * %$\id{from} = \id{to} = 0$% without accessing memory.
   */
  assert(to >= from); if (to == from) return (ALLCLEAR);

  /* This is distressingly complicated.  Here's how it's going to work.
   *
   * There's at least one bit to check, or we'd have returned already --
   * specifically, we must check the bit with index @from@.  But that's at
   * the wrong end.  Instead, we start at the most significant end, with the
   * word containing the bit one short of the @to@ position.  Even though
   * it's actually one off, because we use a half-open interval, we'll call
   * that the `@to@ bit'.
   *
   * We start by loading the word containing the @to@ bit, and start @base@
   * off as the bit index of the least significant bit of this word.  We mask
   * off the high bits (if any), leaving only the @to@ bit and the less
   * significant ones.  We %%\emph{don't}%% check the remaining bits yet.
   *
   * We then start an offset loop.  In each iteration, we check the word
   * we're currently holding: if it's not zero, then we return @base@ plus
   * the position of the most-significant set bit, using @clz32@.  Otherwise,
   * we load the next (less significant) word, and drop @base@ by 32, but
   * don't check it yet.  We end this loop when the word we're holding
   * contains the @from@ bit.  It's possible that we didn't do any iterations
   * of the loop, in which case we're still holding the word containing the
   * @to@ bit at this point.
   *
   * Finally, we mask off the bits below the @from@ bit, and check that what
   * we have left is zero.  If it isn't, we return @base@ plus the position
   * of the most significant bit; if it is, we return @ALLCEAR@.
   */

  /* The first job is to find the word containing the @to@ bit and mask off
   * any higher bits that we don't care about.
   *
   * Recall that the bit's index is @to - 1@, but this must be a valid index
   * because there is at least one bit in the region.  But we start out
   * pointing beyond the vector, so we must add an extra 32 bits.
   */
  n0 = (to + 31)/32; x -= n0; base = (to - 1)&~31u; w = *x++;
  b = to%32; if (b) w &= M32(b);

  /* Before we start the loop, it'd be useful to know how many iterations we
   * need.  This is going to be the offset from the word containing the @to@
   * bit to the word containing the @from@ bit.  Again, we're off by one
   * because that's how our initial indexing is set up.
   */
  n = n0 - from/32 - 1;

  /* Now it's time to do the loop.  This is the easy bit. */
  while (n--) {
    if (w) return (base + 31 - clz32(w));
    w = *x++&MASK32; base -= 32;
  }

  /* We're now holding the final word -- the one containing the @from@ bit.
   * We need to mask off any low bits that we don't care about.
   */
  m = M32(from%32); w &= MASK32&~m;

  /* Final check. */
  if (w) return (base + 31 - clz32(w));
  else return (ALLCLEAR);
}

/*----- General floating-point hacking ------------------------------------*/

/* --- @fltfmt_initbits@ --- *
 *
 * Arguments:	@struct floatbits *x@ = pointer to structure to initialize
 *
 * Returns:	---
 *
 * Use:		Dynamically initialize @x@ to (positive) zero so that it can
 *		be used as the destination operand by other operations.  This
 *		doesn't allocate resources and cannot fail.  The
 *		@FLOATBITS_INIT@ macro is a suitable static initializer for
 *		performing the same task.
 */

void fltfmt_initbits(struct floatbits *x)
{
  x->f = FLTF_ZERO;
  x->a = arena_global;
  x->frac = 0; x->n = x->fracsz = 0;
}

/* --- @fltfmt_freebits@ --- *
 *
 * Arguments:	@struct floatbits *x@ = pointer to structure to free
 *
 * Returns:	---
 *
 * Use:		Releases the memory held by @x@.  Afterwards, @x@ is a valid
 *		(positive) zero, but can safely be discarded.
 */

void fltfmt_freebits(struct floatbits *x)
{
  if (x->frac) x_free(x->a, x->frac);
  x->f = FLTF_ZERO;
  x->frac = 0; x->n = x->fracsz = 0;
}

/* --- @fltfmt_allocfrac@ --- *
 *
 * Arguments:	@struct floatbits *x@ = structure to adjust
 *		@unsigned n@ = number of words required
 *
 * Returns:	---
 *
 * Use:		Reallocate the @frac@ vector so that it has space for at
 *		least @n@ 32-bit words, and set @x->n@ equal to @n@.  If the
 *		current size is already @n@ or greater, then just update the
 *		active length @n@ and return; otherwise, any existing vector
 *		is discarded and a fresh, larger one allocated.
 */

void fltfmt_allocfrac(struct floatbits *x, unsigned n)
  { GROWBUF_REPLACEV(unsigned, x->a, x->frac, x->fracsz, n, 4); x->n = n; }

/* --- @fltfmt_copybits@ --- *
 *
 * Arguments:	@struct floatbits *z_out@ = where to leave the result
 *		@const struct floatbits *x@ = source to copy
 *
 * Returns:	---
 *
 * Use:		Make @z_out@ be a copy of @x@.  If @z_out@ is the same object
 *		as @x@ then do nothing.
 */

void fltfmt_copybits(struct floatbits *z_out, const struct floatbits *x)
{
  unsigned i;

  if (z_out == x) return;
  z_out->f = x->f;
  if (!FINITEP(x)) z_out->exp = 0;
  else z_out->exp = x->exp;
  if ((x->f&(FLTF_ZERO | FLTF_INF)) || !x->n)
    { z_out->n = 0; z_out->frac = 0; }
  else {
    fltfmt_allocfrac(z_out, x->n);
    for (i = 0; i < x->n; i++) z_out->frac[i] = x->frac[i];
  }
}

/* --- @fltfmt_round@ --- *
 *
 * Arguments:	@struct floatbits *z_out@ = destination (may equal source)
 *		@const struct floatbits *x@ = source
 *		@unsigned r@ = rounding mode (@FLTRND_...@ code)
 *		@unsigned n@ = nonzero number of bits to leave
 *
 * Returns:	A @FLTERR_...@ code, specifically either @FLTERR_INEXACT@ if
 *		rounding discarded some nonzero value bits, or @FLTERR_OK@ if
 *		rounding was unnecessary.
 *
 * Use:		Rounds a floating-point value to a given number of
 *		significant bits, using the given rounding rule.
 */

unsigned fltfmt_round(struct floatbits *z_out, const struct floatbits *x,
		      unsigned r, unsigned n)
{
  unsigned rf, i, uw, ub, hw, hb, rc = 0;
  uint32 um, hm, w;
  int exp;

  /* Check that this is likely to work.  We must have at least one bit
   * remaining, so that we can inspect the last-place unit bit.  And we
   * mustn't round up if the current value is already exact, because that
   * would be nonsensical (and inconvenient).
   */
  assert(n > 0); assert(!(r&~(FRPMASK_LOW | FRPMASK_HALF)));

  /* Eliminate trivial cases.  There's nothing to do if the value is infinite
   * or zero, or if we don't have enough precision already.
   *
   * The caller will have set the rounding mode and length suitably for a
   * NaN.
   */
  if (x->f&(FLTF_ZERO | FLTF_INF) || n >= 32*x->n)
    { fltfmt_copybits(z_out, x); return (FLTERR_OK); }

  /* Determine various indices.
   *
   * The quantities @uw@ and @ub@ are the word and bit number which will hold
   * the unit bit when we've finished; @hw@ and @hb@ similarly index the
   * `half' bit, which is the next less significant bit.
   */
	     uw = (n - 1)/32; ub =  -n&31;
  if (!ub) { hw = uw + 1;     hb =     31; }
  else     { hw = uw;	      hb = ub - 1; }

  /* Determine the necessary predicates for the rounding decision. */
  rf = 0;
				  if (x->f&FLTF_NEG) rf |= FRPF_NEG;
  um = B32(ub);			  if (x->frac[uw]&um) rf |= FRPF_ODD;
  hm = B32(hb);			  if (x->frac[hw]&hm) rf |= FRPF_HALF;
				  if (x->frac[hw]&(hm - 1)) rf |= FRPF_LOW;
  for (i = hw + 1; i < x->n; i++) if (x->frac[i]) rf |= FRPF_LOW;
  if (rf&(FRPF_LOW | FRPF_HALF)) rc |= FLTERR_INEXACT;

  /* Allocate space for the result. */
  fltfmt_allocfrac(z_out, uw + 1);

  /* We start looking at the least significant word of the result.  Clear the
   * low bits.
   */
  i = uw; exp = x->exp; w = x->frac[i]&~(um - 1);

  /* If the rounding function is true then we need to add one to the
   * truncated fraction and propagate carries.
   */
  if ((r >> rf)&1) {
    w = (w + um)&MASK32;
    while (i && !w) {
      z_out->frac[i] = w;
      w = (x->frac[--i] + 1)&MASK32;
    }
    if (!w) { w = B31; exp++; }
  }

  /* Store, and copy the remaining words. */
  for (;;) {
    z_out->frac[i] = w;
    if (!i) break;
    w = x->frac[--i];
  }

  /* Done. */
  z_out->f = x->f&(FLTF_NEG | FLTF_NANMASK);
  if (x->f&FLTF_NANMASK) z_out->exp = 0;
  else z_out->exp = exp;
  return (rc);
}

/*----- IEEE and related formats ------------------------------------------*/

/* IEEE (and related) format descriptions. */
const struct fltfmt_ieeefmt
  fltfmt_mini = { FLTIF_HIDDEN, 4, 4 },
  fltfmt_bf16 = { FLTIF_HIDDEN, 8, 8 },
  fltfmt_f16 = { FLTIF_HIDDEN, 5, 11 },
  fltfmt_f32 = { FLTIF_HIDDEN, 8, 24 },
  fltfmt_f64 = { FLTIF_HIDDEN, 11, 53 },
  fltfmt_f128 = { FLTIF_HIDDEN, 15, 113 },
  fltfmt_idblext80 = { 0, 15, 64 };

/* --- @fltfmt_encieee@ ---
 *
 * Arguments:	@const struct fltfmt_ieeefmt *fmt@ = format description
 *		@uint32 *z@ = output vector
 *		@const struct floatbits *x@ = value to encode
 *		@unsigned r@ = rounding mode
 *		@unsigned errmask@ = error mask
 *
 * Returns:	Error flags (@FLTERR_...@).
 *
 * Use:		Encode a floating-point value in an IEEE format.  This is the
 *		machinery shared by the @fltfmt_enc...@ functions for
 *		encoding IEEE-format values.  Most of the arguments and
 *		behaviour are as described for those functions.
 *
 *		The encoded value is right-aligned and big-endian; i.e., the
 *		sign bit ends up in @z[0]@, and the least significant bit of
 *		the significand ends up in the least significant bit of
 *		@z[n - 1]@.
 */

unsigned fltfmt_encieee(const struct fltfmt_ieeefmt *fmt,
			uint32 *z, const struct floatbits *x,
			unsigned r, unsigned errmask)
{
  struct floatbits y;
  unsigned sigwd, fracwd, err = 0, f = x->f, rf;
  unsigned i, j, n, nb, nw, mb, mw, esh, sh;
  int exp, minexp, maxexp;
  uint32 z0, t;

#define ERR(e) do { err |= (e); if (err&~errmask) goto end; } while (0)

  /* The following code assumes that the sign, biased exponent, unit, and
   * quiet/signalling bits can all fit into the most significant 32 bits of
   * the result.
   */
  assert(fmt->expwd + 3 <= 32);
  esh = 31 - fmt->expwd;

  /* Determine the output size. */
  nb = fmt->prec + fmt->expwd + 1;
  if (fmt->f&FLTIF_HIDDEN) nb--;
  nw = (nb + 31)/32;

  /* Determine the top bits. */
  z0 = 0; i = 0;
  if (x->f&FLTF_NEG) z0 |= B31;

  /* And now for the main case analysis. */

  if (f&FLTF_ZERO) {
    /* Zero.  There's very little to do here. */

  } else if (f&FLTF_INF) {
    /* Infinity.  Set the exponent and, for non-hidden-bit formats, the unit
     * bit.
     */

    z0 |= M32(fmt->expwd) << esh;
    if (!(fmt->f&FLTIF_HIDDEN)) z0 |= B32(esh - 1);

  } else if (f&FLTF_NANMASK) {
    /* Not-a-number.
     *
     * We must check that we won't lose significant bits.  We need a bit for
     * the quiet/signalling flag, and enough space for the significant
     * payload bits.  The unit bit is not in play here, so the available
     * space is always one less than the advertised precision.  To put it
     * another way, we need space for the payload, a bit for the
     * quiet/signalling flag, and a bit for the unit place.
     */

    fracwd = sigbits(x);
    if (fracwd + 2 > fmt->prec) ERR(FLTERR_INEXACT);

    /* Copy the payload.
     *
     * If the payload is all-zero and we're meant to set a signalling NaN
     * then report an exactness failure and set the least-significant bit.
     */
    mb = fmt->prec - 2; mw = (mb + 31)/32; sh = -mb%32;
    n = x->n;
      if (n < mw) j = 0;
      else { n = mw; j = sh; }
    if ((f&FLTF_SNAN) && ms_set_bit(x->frac + n, j, 32*n) == ALLCLEAR) {
      ERR(FLTERR_INEXACT);
      n = nw - 1; for (i = 0; i < n; i++) z[i] = 0;
      z[i++] = 1;
    } else {
      for (i = 0; i < nw - mw; i++) z[i] = 0;
      n = x->n; if (n > mw) n = mw;
      t = shr(z + i, x->frac, n, sh); i += n;
      if (i < nw) z[i++] = t;
      sh = esh - 2; if (fmt->f&FLTIF_HIDDEN) sh++;
      if (f&FLTF_QNAN) z0 |= B32(sh);
    }

    /* Set the exponent and, for non-hidden-bit formats, the unit bit. */
    z0 |= M32(fmt->expwd) << esh;
    if (!(fmt->f&FLTIF_HIDDEN)) z0 |= B32(esh - 1);

  } else {
    /* A finite value.
     *
     * Adjust the exponent by one place to compensate for the difference in
     * significant conventions.  Our significand lies between zero (in fact,
     * a half, because we require normalization) and one, while an IEEE
     * significand lies between zero (in fact, one) and two.  Our exponent is
     * therefore one larger than the IEEE exponent will be.
     */

    /* Determine the maximum true (unbiased) exponent.  As noted above, this
     * is also the bias.
     */
    exp = x->exp - 1;
    maxexp = (1 << (fmt->expwd - 1)) - 1;
    minexp = 1 - maxexp;

    if (exp <= minexp - (int)fmt->prec) {
      /* If the exponent is very small then we underflow.  We have %$p - 1$%
       * bits available to represent a subnormal significand, and therefore
       * can represent at least one bit of a value as small as
       * %$2^{e_{\text{min}}-p+1}$%.
       *
       * If the exponent is one short of the threshold, then we check to see
       * whether the value will round up.
       */

      if ((minexp - exp == fmt->prec) &&
	  ((r >> (FRPF_HALF |
		  (sigbits(x) > 1 ? FRPF_LOW : 0) |
		  (f&FLTF_NEG ? FRPF_NEG : 0)))&1)) {
	ERR(FLTERR_INEXACT);
	for (i = 0; i < nw - 1; i++) z[i] = 0;
	z[i++] = 1;
      } else {
	ERR(FLTERR_UFLOW | FLTERR_INEXACT);
	/* Return (signed) zero. */
      }

    } else {
      /* We can at least try to store some bits. */

      /* Let's see how many we need to deal with and how much space we have.
       * We might as well set the biased exponent here while we're at it.
       *
       * If %$e \ge e_{\text{min}}$% then we can store %$p$% bits of
       * significand.  Otherwise, we must make a subnormal and we can only
       * store %$p + e - e_{\text{min}}$% bits.  (Cross-check: if %$e \le
       * e_{\text{min}} - p$% then we can store zero bits or fewer and have
       * underflowed to zero, which matches the previous case.)  In the
       * subnormal case, we also `correct' the exponent so that we store the
       * correct sentinel value later.
       */
      fracwd = sigbits(x);
      if (exp >= minexp) sigwd = fmt->prec;
      else { sigwd = fmt->prec + exp - minexp; exp = minexp - 1; }
      mw = (sigwd + 31)/32; sh = -sigwd%32;

      /* If we don't have enough significand bits then we must round.  This
       * might increase the exponent, so we must reload.
       */
      if (fracwd > sigwd) {
	ERR(FLTERR_INEXACT);
	y.frac = z + nw - mw; y.fracsz = mw; fltfmt_round(&y, x, r, sigwd);
	x = &y; exp = y.exp - 1; fracwd = sigwd;
      }

      if (exp > maxexp) {
	/* If the exponent is too large, then we overflow.  If the error is
	 * masked, then we must produce a default value, choosing between
	 * infinity and the largest representable finite value according to
	 * the rounding mode.
	 */

	ERR(FLTERR_OFLOW | FLTERR_INEXACT);
	rf = FRPF_ODD | FRPF_HALF | FRPF_LOW;
	if (f&FLTF_NEG) rf |= FRPF_NEG;
	if ((r >> rf)&1)
	  z0 |= M32(fmt->expwd) << esh;
	else {
	  z0 |= (B32(fmt->expwd) - 2) << esh;
	  mb = fmt->prec; if (fmt->f&FLTIF_HIDDEN) mb--;
	  mw = (mb + 31)/32;
	  i = nw - mw;
	  z[i++] = M32(mb%32);
	  while (i < nw) z[i++] = MASK32;
	}

      } else {
	/* The exponent is in range.  Everything is ready. */

	/* Store the significand. */
	n = (fracwd + 31)/32; i = nw - mw;
	t = shr(z + i, x->frac, n, sh); i += n;
	if (i < nw) z[i++] = t;

	/* Fill in the top end. */
	for (j = nw - mw; j--; ) z[j] = 0;

	/* Set the biased exponent. */
	z0 |= (exp + maxexp) << esh;

	/* Clear the unit bit if we're suppose to use a hidden-bit convention. */
	if (fmt->f&FLTIF_HIDDEN) {
	  mb = fmt->prec - 1; mw = (mb + 31)/32; mb = mb%32;
	  z[nw - mw] &= ~B32(mb);
	}
      }
    }
  }

  /* Clear the significand bits that we haven't set explicitly yet. */
  while (i < nw) z[i++] = 0;

  /* All that remains is to insert the top bits @z0@ in the right place.
   * This will set the exponent, and the unit and quiet bits.
   */
  sh = -nb%32;
  z[0] |= z0 >> sh;
  if (sh && nb >= 32) z[1] |= z0 << (32 - sh);

end:
  return (err);

#undef ERR
}

/* --- @fltfmt_encTY@ --- *
 *
 * Arguments:	@octet *z_out@, @uint16 *z_out@, @uint32 *z_out@,
 *			@kludge64 *z_out@ = where to put the encoded value
 *		@uint16 *se_out@, @kludge64 *m_out@ = where to put the
 *			encoded sign-and-exponent and significand
 *		@const struct floatbits *x@ = value to encode
 *		@unsigned r@ = rounding mode
 *		@unsigned errmask@ = error mask
 *
 * Returns:	Error flags (@FLTERR_...@).
 *
 * Use:		Encode a floating-point value in an IEEE (or IEEE-adjacent)
 *		format.
 *
 *		If an error is encountered during the encoding, and the
 *		corresponding bit of @errmask@ is clear, then processing
 *		stops immediately and the error is returned; if the bit is
 *		set, then processing continues as described below.
 *
 *		The @TY@ may be
 *
 *		  * @mini@ for the 8-bit `1.4.3 minifloat' format, with
 *		    four-bit exponent and four-bit significand, represented
 *		    as a single octet;
 *
 *		  * @bf16@ for the Google `bfloat16' format, with eight-bit
 *		    exponent and eight-bit significand, represented as a
 *		    @uint16@;
 *
 *		  * @f16@ for the IEEE `binary16' format, with five-bit
 *		    exponent and eleven-bit significand, represented as a
 *		    @uint16@;
 *
 *		  * @f32@ for the IEEE `binary32' format, with eight-bit
 *		    exponent and 24-bit significand, represented as a
 *		    @uint32@;
 *
 *		  * @f64@ for the IEEE `binary64' format, with eleven-bit
 *		    exponent and 53-bit significand, represented as a
 *		    @kludge64@;
 *
 *		  * @f128@ for the IEEE `binary128' format, with fifteen-bit
 *		    exponent and 113-bit significand, represented as four
 *		    @uint32@ limbs, most significant first; or
 *
 *		  * @idblext80@ for the Intel 80-bit `double extended'
 *		    format, with fifteen-bit exponent and 64-bit significand
 *		    with no hidden bit, represented as a @uint16 se@
 *		    holding the sign and exponent, and a @kludge64 m@
 *		    holding the significand.
 *
 *		Positive and negative zero and infinity are representable
 *		exactly.
 *
 *              Following IEEE recommendations (and most implementations),
 *		the most significant fraction bit of a quiet NaN is set; this
 *		bit is clear in a signalling NaN.  The most significant
 *		payload bits of a NaN, held in the top bits of @x->frac[0]@,
 *		are encoded in the output significand following the `quiet'
 *		bit.  If the chosen format's significand field is too small
 *		to accommodate all of the set payload bits then the
 *		@FLTERR_INEXACT@ error bit is set and, if masked, the
 *		excess payload bits are discarded.  No rounding of NaN
 *		payloads is performed.
 *
 *		Otherwise, the input value is finite and nonzero.  If the
 *		significand cannot be represented exactly then the
 *		@FLTERR_INEXACT@ error bit is set, and, if masked, the value
 *		will be rounded (internally -- the input @x@ is not changed).
 *		If the (rounded) value's exponent is too large to represent,
 *		then the @FLTERR_OFLOW@ and @FLTERR_INEXACT@ error bits are
 *		set and, if masked, the result is either the (absolute)
 *		largest representable finite value or infinity, with the
 *		appropriate sign, chosen according to the rounding mode.  If
 *		the exponent is too small to represent, then the
 *		@FLTERR_UFLOW@ and @FLTERR_INEXACT@ error bits are set and,
 *		if masked, the result is either the (absolute) smallest
 *		nonzero value or zero, with the appropriate sign, chosen
 *		according to the rounding mode.
 */

unsigned fltfmt_encmini(octet *z_out, const struct floatbits *x,
			unsigned r, unsigned errmask)
{
  uint32 t[1];
  unsigned rc;

  rc = fltfmt_encieee(&fltfmt_mini, t, x, r, errmask);
  if (!(rc&~errmask)) *z_out = t[0];
  return (rc);
}

unsigned fltfmt_encbf16(uint16 *z_out, const struct floatbits *x,
			unsigned r, unsigned errmask)
{
  uint32 t[1];
  unsigned rc;

  rc = fltfmt_encieee(&fltfmt_bf16, t, x, r, errmask);
  if (!(rc&~errmask)) *z_out = t[0];
  return (rc);
}

unsigned fltfmt_encf16(uint16 *z_out, const struct floatbits *x,
		       unsigned r, unsigned errmask)
{
  uint32 t[1];
  unsigned rc;

  rc = fltfmt_encieee(&fltfmt_f16, t, x, r, errmask);
  if (!(rc&~errmask)) *z_out = t[0];
  return (rc);
}

unsigned fltfmt_encf32(uint32 *z_out, const struct floatbits *x,
		       unsigned r, unsigned errmask)
  { return (fltfmt_encieee(&fltfmt_f32, z_out, x, r, errmask)); }

unsigned fltfmt_encf64(kludge64 *z_out, const struct floatbits *x,
		       unsigned r, unsigned errmask)
{
  uint32 t[2];
  unsigned rc;

  rc = fltfmt_encieee(&fltfmt_f64, t, x, r, errmask);
  if (!(rc&~errmask)) SET64(*z_out, t[0], t[1]);
  return (rc);
}

unsigned fltfmt_encf128(uint32 *z_out, const struct floatbits *x,
			unsigned r, unsigned errmask)
  { return (fltfmt_encieee(&fltfmt_f128, z_out, x, r, errmask)); }

unsigned fltfmt_encidblext80(uint16 *se_out, kludge64 *m_out,
			     const struct floatbits *x,
			     unsigned r, unsigned errmask)
{
  uint32 t[3];
  unsigned rc;

  rc = fltfmt_encieee(&fltfmt_idblext80, t, x, r, errmask);
  if (!(rc&~errmask)) { *se_out = t[0]; SET64(*m_out, t[1], t[2]); }
  return (rc);
}

/* --- @fltfmt_decieee@ --- *
 *
 * Arguments:	@const struct fltfmt_ieeefmt *fmt@ = format description
 *		@struct floatbits *z_out@ = output decoded representation
 *		@const uint32 *x@ = input encoding
 *
 * Returns:	Error flags (@FLTERR_...@).
 *
 * Use:		Decode a floating-point value in an IEEE format.  This is the
 *		machinery shared by the @fltfmt_dec...@ functions for
 *		deccoding IEEE-format values.  Most of the arguments and
 *		behaviour are as described for those functions.
 *
 *		The encoded value should be right-aligned and big-endian;
 *		i.e., the sign bit ends up in @z[0]@, and the least
 *		significant bit of the significand ends up in the least
 *		significant bit of @z[n - 1]@.
 */

unsigned fltfmt_decieee(const struct fltfmt_ieeefmt *fmt,
			struct floatbits *z_out, const uint32 *x)
{
  unsigned sigwd, err = 0, f = 0;
  unsigned i, nb, nw, mw, esh, sh;
  int exp, minexp, maxexp;
  uint32 x0, t, u, emask;

  /* The following code assumes that the sign, biased exponent, unit, and
   * quiet/signalling bits can all fit into the most significant 32 bits of
   * the result.
   */
  assert(fmt->expwd + 3 <= 32);
  esh = 31 - fmt->expwd; emask = M32(fmt->expwd);
  sigwd = fmt->prec; if (fmt->f&FLTIF_HIDDEN) sigwd--;

  /* Determine the input size. */
  nb = sigwd + fmt->expwd + 1; nw = (nb + 31)/32;

  /* Extract the sign, exponent, and top of the significand. */
  sh = -nb%32;
  x0 = x[0] << sh;
  if (sh && nb >= 32) x0 |= x[1] >> (32 - sh);
  if (x0&B31) f |= FLTF_NEG;
  t = (x0 >> esh)&emask;

  /* Time for a case analysis. */

  if (t == emask) {
    /* Infinity or NaN.
     *
     * Note that we don't include the quiet bit in our decoded payload.
     */

    if (!(fmt->f&FLTIF_HIDDEN)) {
      /* No hidden bit, so we expect the unit bit to be set.  If it isn't,
       * that's technically invalid, and its absence won't survive a round
       * trip, since the bit isn't considered part of a NaN payload -- or
       * even to distinguish a NaN from an infinity.  In any event, reduce
       * the notional significand size to exclude this bit from further
       * consideration.
       */

      if (!(x0&B32(esh - 1))) err = FLTERR_INVAL;
      sigwd--;
    }

    if (ms_set_bit(x + nw, 0, sigwd) == ALLCLEAR)
      f |= FLTF_INF;
    else {
      sh = esh - 2; if (fmt->f&FLTIF_HIDDEN) sh++;
      if (x0&B32(sh)) f |= FLTF_QNAN;
      else f |= FLTF_SNAN;
      sigwd--; mw = (sigwd + 31)/32;
      fltfmt_allocfrac(z_out, mw);
      shl(z_out->frac, x + nw - mw, mw, -sigwd%32);
    }
    goto end;
  }

  /* Determine the exponent bounds. */
  maxexp = (1 << (fmt->expwd - 1)) - 1;
  minexp = 1 - maxexp;

  /* Dispatch.  If there's a hidden bit then everything is well defined.
   * Otherwise, we'll normalize the incoming value regardless, but report
   * settings of the unit bit which are inconsistent with the exponent.
   */
  if (fmt->f&FLTIF_HIDDEN) {
    if (!t) { exp = minexp; goto normalize; }
    else { exp = t - maxexp; goto hidden; }
  } else {
    u = x0&B32(esh - 1);
    if (!t) { exp = minexp; if (u) err |= FLTERR_INVAL; }
    else { exp = t - maxexp; if (!u) err |= FLTERR_INVAL; }
    goto normalize;
  }

hidden:
  /* We have a normal real number with a hidden bit. */

  mw = (sigwd + 31)/32;

  if (!(sigwd%32)) {
    /* The bits we have occupy a whole number of words, but we need to shift
     * to make space for the unit bit.
     */

    fltfmt_allocfrac(z_out, mw + 1);
    z_out->frac[mw] = shr(z_out->frac, x + nw - mw, mw, 1);
  } else {
    fltfmt_allocfrac(z_out, mw);
    shl(z_out->frac, x + nw - mw, mw, -(sigwd + 1)%32);
  }
  z_out->frac[0] |= B31;
  z_out->exp = exp + 1;
  goto end;

normalize:
  /* We have, at least potentially, a subnormal number, with no hidden
   * bit.
   */

  i = ms_set_bit(x + nw, 0, sigwd);
  if (i == ALLCLEAR) { f |= FLTF_ZERO; goto end; }
  mw = i/32 + 1; sh = 32*mw - i - 1;
  fltfmt_allocfrac(z_out, mw);
  shl(z_out->frac, x + nw - mw, mw, sh);
  z_out->exp = exp - fmt->prec + 2 + i;
  goto end;

end:
  /* All done. */
  z_out->f = f; return (err);
}

/* --- @fltfmt_decTY@ --- *
 *
 * Arguments:	@const struct floatbits *z_out@ = storage for the result
 *		@octet x@, @uint16 x@, @uint32 x@, @kludge64 x@ =
 *			encoded input
 *		@uint16 se@, @kludge64 m@ = encoded sign-and-exponent and
 *			significand
 *
 * Returns:	Error flags (@FLTERR_...@).
 *
 * Use:		Encode a floating-point value in an IEEE (or IEEE-adjacent)
 *		format.
 *
 *		The options for @TY@ are as documented for the encoding
 *		functions above.
 *
 *		In formats without a hidden bit -- currently only @idblext80@
 *		-- not all bit patterns are valid encodings.  If the explicit
 *		unit bit is set when the exponent field is all-bits-zero, or
 *		clear when the exponent field is not all-bits-zero, then the
 *		@FLTERR_INVAL@ error bit is set.  If the exponent is all-
 *		bits-set, denoting infinity or a NaN, then the unit bit is
 *		otherwise ignored -- in particular, it does not affect the
 *		NaN payload, or even whether the input encodes a NaN or
 *		infinity.  Otherwise, the unit bit is considered significant,
 *		and the result is normalized as one would expect.
 *		Consequently, biased exponent values 0 and 1 are distinct
 *		only with respect to which bit patterns are considered valid,
 *		and not with respect to the set of values denoted.
 */

unsigned fltfmt_decmini(struct floatbits *z_out, octet x)
  { uint32 t[1]; t[0] = x; return (fltfmt_decieee(&fltfmt_mini, z_out, t)); }

unsigned fltfmt_decbf16(struct floatbits *z_out, uint16 x)
  { uint32 t[1]; t[0] = x; return (fltfmt_decieee(&fltfmt_bf16, z_out, t)); }

unsigned fltfmt_decf16(struct floatbits *z_out, uint16 x)
  { uint32 t[1]; t[0] = x; return (fltfmt_decieee(&fltfmt_f16, z_out, t)); }

unsigned fltfmt_decf32(struct floatbits *z_out, uint32 x)
  { uint32 t[1]; t[0] = x; return (fltfmt_decieee(&fltfmt_f32, z_out, t)); }

unsigned fltfmt_decf64(struct floatbits *z_out, kludge64 x)
{
  uint32 t[2];

  t[0] = HI64(x); t[1] = LO64(x);
  return (fltfmt_decieee(&fltfmt_f64, z_out, t));
}

unsigned fltfmt_decf128(struct floatbits *z_out, const uint32 *x)
  { return (fltfmt_decieee(&fltfmt_f128, z_out, x)); }

unsigned fltfmt_decidblext80(struct floatbits *z_out, uint16 se, kludge64 m)
{
  uint32 t[3];

  t[0] = se; t[1] = HI64(m); t[2] = LO64(m);
  return (fltfmt_decieee(&fltfmt_idblext80, z_out, t));
}

/*----- Native formats ----------------------------------------------------*/

/* If the floating-point radix is a power of two, determine how many bits
 * there are in each digit.  This isn't exhaustive, but it covers most of the
 * bases, so to speak.
 */
#if FLT_RADIX == 2
#  define DIGIT_BITS 1
#elif FLT_RADIX == 4
#  define DIGIT_BITS 2
#elif FLT_RADIX == 8
#  define DIGIT_BITS 3
#elif FLT_RADIX == 16
#  define DIGIT_BITS 4
#endif

/* Take note if we need to cope with the revered quiet/signalling convention
 * used by HP-PA and older MIPS processors.
 */
#if defined(__hppa__) || (defined(__mips__) && !defined(__mips_nan2008))
#  define FROB_NANS
#endif

/* --- @ENCFLT@ --- *
 *
 * Arguments:	@ty@ = the C type to encode
 *		@TY@ = the uppercase prefix for macros
 *		@ty (*ldexp)(ty, int)@ = function to scale a @ty@ value by a
 *			power of two
 *		@unsigned &rc@ = error code to set
 *		@ty *z_out@ = storage for the result
 *		@const struct floatbits *x@ = value to convert
 *		@unsigned r@ = rounding mode
 *
 * Returns:	---
 *
 * Use:		Encode a floating-point value @x@ as a native C object of
 *		type @ty@.  This is the machinery shared by the
 *		@fltfmt_enc...@ functions for enccoding native-format values.
 *		Most of the behaviour is as described for those functions.
 */

/* Utilities based on conditional compilation, because we can't use
 * %|#ifdef|% directly in macros.
 */

#ifdef NAN
   /* The C implementation acknowledges the existence of (quiet) NaN values,
    * but will neither let us set the payload in a useful way, nor
    * acknowledge the existence of signalling NaNs.  We have no good way to
    * determine which NaN the @NAN@ macro produces, so report this conversion
    * as inexact.
    */

#  define SETNAN(rc, z) do { (z) = NAN; (rc) = FLTERR_INEXACT; } while (0)
#else
   /* This C implementation doesn't recognize NaNs.  This value is totally
    * unrepresentable, so just report the error.  (Maybe it's C89 and would
    * actually do the right thing with @0/0@.  I'm not sure the requisite
    * compile-time configuration machinery is worth the effort.)
    */

#  define SETNAN(rc, z) do { (z) = 0; (rc) = FLTERR_REPR; } while (0)
#endif

#ifdef INFINITY
   /* The C implementation supports infinities.  This is a simple win. */

#  define SETINF(TY, rc, z)						\
	do { (z) = INFINITY; (rc) = FLTERR_OK; } while (0)
#else
   /* The C implementation doesn't support infinities.  Return the maximum
    * value and report it as an overflow; I think this is more useful than
    * reporting a complete representation failure.  (Maybe it's C89 and would
    * actually do the right thing with @1/0@.  Again, I'm not sure the
    * requisite compile-time configuration machinery is worth the effort.)
    */

#  define SETINF(TY, rc, z)						\
	do { (z) = TY##_MAX; (rc) = FLTERR_OFLOW | FLTERR_INEXACT; } while (0)
#endif

#ifdef DIGIT_BITS
   /* The floating point formats use a power-of-two radix.  This means that
    * we can determine the correctly rounded value before we start building
    * the native floating-point value.
    */

#  define ENC_ROUND_DECLS struct floatbits _y;
#  define ENC_ROUND(TY, rc, x, r) do {					\
     (rc) |= fltfmt_round(&_y, (x), (r), DIGIT_BITS*TY##_MANT_DIG);	\
     (x) = &_y;								\
   } while (0)
#else
   /* The floating point formats use a non-power-of-two radix.  This means
    * that conversion is inherently inexact.
    */

#  define ENC_ROUND_DECLS
#  define ENC_ROUND(TY, rc, x, r)					\
     do (rc) |= FLTERR_INEXACT; while (0)
#  define ENC_FIXUP(...)

#endif

#ifdef FROB_NANS
#  define FROBNAN_ENCDECLS	struct floatbits _y
#  define FROBNAN_ENC do {						\
     if (_x->f&FLTF_NANMASK) {						\
       _y.f = _x->f ^ FLTF_NANMASK; _y.frac = _x->frac; _y.n = _x->n;	\
       _x = &_y;							\
     }									\
   } while (0)
#else
#  define FROBNAN_ENCDECLS
#  define FROBNAN_ENC do ; while (0)
#endif

#define ENCFLT(ty, TY, ldexp, rc, z_out, x, r) do {			\
  const struct floatbits *_x = (x);					\
  unsigned _rc = 0;							\
  FROBNAN_ENCDECLS;							\
									\
  /* See if the native format is one that we recognize. */		\
  switch (TY##_FORMAT&(FLTFMT_ORGMASK | FLTFMT_TYPEMASK)) {		\
									\
    case FLTFMT_IEEE_F32: {						\
      uint32 _t[1];							\
      unsigned char *_z = (unsigned char *)(z_out);			\
									\
      FROBNAN_ENC;							\
      (rc) = fltfmt_encieee(&fltfmt_f32, _t, _x, (r), FLTERR_ALLERRS);	\
      switch (TY##_FORMAT&FLTFMT_ENDMASK) {				\
	case FLTFMT_BE: STORE32_B(_z, _t[0]); break;			\
	case FLTFMT_LE: STORE32_L(_z, _t[0]); break;			\
	default: assert(!"unimplemented byte order"); break;		\
      }									\
    } break;								\
									\
    case FLTFMT_IEEE_F64: {						\
      uint32 _t[2];							\
      unsigned char *_z = (unsigned char *)(z_out);			\
									\
      FROBNAN_ENC;							\
      (rc) = fltfmt_encieee(&fltfmt_f64, _t, _x, (r), FLTERR_ALLERRS);	\
      switch (TY##_FORMAT&FLTFMT_ENDMASK) {				\
	case FLTFMT_BE:							\
	  STORE32_B(_z + 0, _t[0]); STORE32_B(_z + 4, _t[1]);		\
	  break;							\
	case FLTFMT_LE:							\
	  STORE32_L(_z + 0, _t[1]); STORE32_L(_z + 4, _t[0]);		\
	  break;							\
	case FLTFMT_ARME:						\
	  STORE32_L(_z + 0, _t[0]); STORE32_L(_z + 4, _t[1]);		\
	  break;							\
	default: assert(!"unimplemented byte order"); break;		\
      }									\
    } break;								\
									\
    case FLTFMT_IEEE_F128: {						\
      uint32 _t[4];							\
      unsigned char *_z = (unsigned char *)(z_out);			\
									\
      FROBNAN_ENC;							\
      (rc) = fltfmt_encieee(&fltfmt_f128, _t, _x, (r), FLTERR_ALLERRS);	\
      switch (TY##_FORMAT&FLTFMT_ENDMASK) {				\
	case FLTFMT_BE:							\
	  STORE32_B(_z +  0, _t[0]); STORE32_B(_z +  4, _t[1]);		\
	  STORE32_B(_z +  8, _t[0]); STORE32_B(_z + 12, _t[1]);		\
	  break;							\
	case FLTFMT_LE:							\
	  STORE32_L(_z +  0, _t[3]); STORE32_L(_z +  4, _t[2]);		\
	  STORE32_L(_z +  8, _t[1]); STORE32_L(_z + 12, _t[0]);		\
	  break;							\
	default: assert(!"unimplemented byte order"); break;		\
      }									\
    } break;								\
									\
    case FLTFMT_INTEL_F80: {						\
      uint32 _t[3];							\
      unsigned char *_z = (unsigned char *)(z_out);			\
									\
      FROBNAN_ENC;							\
      (rc) = fltfmt_encieee(&fltfmt_idblext80,				\
			    _t, _x, (r), FLTERR_ALLERRS);		\
      switch (TY##_FORMAT&FLTFMT_ENDMASK) {				\
	case FLTFMT_BE:							\
	  STORE16_B(_z + 0, _t[0]);					\
	  STORE32_B(_z + 2, _t[1]); STORE32_B(_z + 6, _t[2]);		\
	  break;							\
	case FLTFMT_LE:							\
	  STORE32_L(_z + 0, _t[2]); STORE32_L(_z + 4, _t[1]);		\
	  STORE16_L(_z + 8, _t[0]);					\
	  break;							\
	default: assert(!"unimplemented byte order"); break;		\
      }									\
    } break;								\
									\
    default: {								\
      /* We must do this the hard way. */				\
									\
      ty _z;								\
      unsigned _i;							\
      ENC_ROUND_DECLS;							\
									\
      /* Case analysis... */						\
      if (_x->f&FLTF_NANMASK) {						\
	/* A NaN.  Use the macro above. */				\
									\
	SETNAN(_rc, _z);						\
	if (x->f&FLTF_NEG) _z = -_z;					\
      } else if (_x->f&FLTF_INF) {					\
	/* Infinity.  Use the macro. */					\
									\
	SETINF(TY, _rc, _z);						\
	if (_x->f&FLTF_NEG) _z = -_z;					\
      } else if (_x->f&FLTF_ZERO) {					\
	/* Zero.  If we're asked for a negative zero then check that we	\
	 * produced one: if not, then report an exactness failure.	\
	 */								\
									\
	_z = 0.0;							\
	if (_x->f&FLTF_NEG)						\
	  { _z = -_z; if (!NEGP(_z)) _rc |= FLTERR_INEXACT; }		\
      } else {								\
	/* A finite value. */						\
									\
	/* If the radix is a power of two, we can round to the correct	\
	 * precision, which will save inexactness later.		\
	 */								\
	ENC_ROUND(TY, _rc, _x, (r));					\
									\
	/* Insert the 32-bit pieces of the fraction one at a time,	\
	 * starting from the least-significant end.  This minimizes the	\
	 * inaccuracy from the overall approach, but it's imperfect	\
	 * unless the value has already been rounded correctly.		\
	 */								\
	_z = 0.0;							\
	for (_i = _x->n, _z = 0.0; _i--; )				\
	  _z += ldexp(_x->frac[_i], _x->exp - 32*_i);			\
									\
	/* Negate the value if we need to. */				\
	if (_x->f&FLTF_NEG) _z = -_z;					\
      }									\
									\
      /* All done. */							\
      *(z_out) = _z;							\
    } break;								\
  }									\
									\
  /* Set the error code. */						\
  (rc) = _rc;								\
} while (0)

/* --- @fltfmt_encTY@ --- *
 *
 * Arguments:	@ty *z_out@ = storage for the result
 *		@const struct floatbits *x@ = value to encode
 *		@unsigned r@ = rounding mode
 *
 * Returns:	Error flags (@FLTERR_...@).
 *
 * Use:		Encode the floating-point value @x@ as a native C object and
 *		store the result in @z_out@.
 *
 *		The @TY@ may be @flt@ to encode a @float@, @dbl@ to encode a
 *		@double@, or (on C99 implementations) @ldbl@ to encode a
 *		@long double@.
 *
 *		In detail, conversion is performed as follows.
 *
 *		  * If a non-finite value cannot be represented by the
 *		    implementation then the @FLTERR_REPR@ error bit is set
 *		    and @*z_out@ is set to zero if @x@ is a NaN, or the
 *		    (absolute) largest representable value, with appropriate
 *		    sign, if @x@ is an infinity.
 *
 *		  * If the implementation can represent NaNs, but cannot set
 *		    NaN payloads, then the @FLTERR_INEXACT@ error bit is set,
 *		    and @*z_out@ is set to an arbitrary (quiet) NaN value.
 *
 *		  * If @x@ is negative zero, but the implementation does not
 *		    distinguish negative and positive zero, then the
 *		    @FLTERR_INEXACT@ error bit is set and @*z_out@ is set to
 *		    zero.
 *
 *		  * If the implementation's floating-point radix is not a
 *		    power of two, and @x@ is a nonzero finite value, then
 *		    @FLTERR_INEXACT@ error bit is set (unconditionally), and
 *		    the value is rounded by the implementation using its
 *		    prevailing rounding policy.  If the radix is a power of
 *		    two, then the @FLTERR_INEXACT@ error bit is set only if
 *		    rounding is necessary, and rounding is performed using
 *		    the rounding mode @r@.
 */

unsigned fltfmt_encflt(float *z_out, const struct floatbits *x, unsigned r)
{
  unsigned rc;

  ENCFLT(double, FLT, ldexp, rc, z_out, x, r);
  return (rc);
}

unsigned fltfmt_encdbl(double *z_out, const struct floatbits *x, unsigned r)
{
  unsigned rc;

  ENCFLT(double, DBL, ldexp, rc, z_out, x, r);
  return (rc);
}

#if __STDC_VERSION__ >= 199001
unsigned fltfmt_encldbl(long double *z_out,
			const struct floatbits *x, unsigned r)
{
  unsigned rc;

  ENCFLT(long double, LDBL, ldexpl, rc, z_out, x, r);
  return (rc);
}
#endif

/* --- @DECFLT@ --- *
 *
 * Arguments:	@ty@ = the C type to encode
 *		@TY@ = the uppercase prefix for macros
 *		@ty (*frexp)(ty, int *)@ = function to decompose a @ty@ value
 *			into a binary exponent and normalized fraction
 *		@unsigned &rc@ = error code to set
 *		@struct floatbits *z_out@ = storage for the result
 *		@ty x@ = value to convert
 *		@unsigned r@ = rounding mode
 *
 * Returns:	---
 *
 * Use:		Decode a C native floating-point object.  This is the
 *		machinery shared by the @fltfmt_dec...@ functions for
 *		decoding native-format values.  Most of the behaviour is as
 *		described for those functions.
 */

/* Define some utilities for decoding native floating-point formats.
 *
 *   * @NFRAC(d)@ is the number of fraction limbs we'll need for @d@ native
 *     digits.
 *
 *   * @CONVFIX@ is a final fixup applied to the decoded value.
 */
#ifdef DIGIT_BITS
#  define NFRAC(TY) ((DIGIT_BITS*TY##_MANT_DIG + 31)/32)
#  define CONVFIX(ty, rc, z, x, n, f, r) do assert(!(x)); while (0)
#else
#  define NFRAC(TY)							\
     (ceil(log(pow(FLT_RADIX, TY##_MANT_DIG) - 1)/32.0*log(2.0)) + 1)
#  define CONVFIX(ty, rc, z, x, n, f, r) do {				\
     ty _x_ = (x);							\
     struct floatbits *_z_ = (z);					\
     uint32 _w_;							\
     unsigned _i_, _n_ = (n), _f_;					\
									\
     /* Round the result according to the remainder left in %$x$%. */	\
     _f_ = 0; _i_ = _n_ - 1; _w_ = _z_->frac[_i_];			\
     if ((f)&FLTF_NEG) _f_ |= FRPF_NEG;					\
     if (_w_&1) _f_ |= FRPF_ODD;					\
     if (_y_ >= 0.5) _f_ |= FRPF_HALF;					\
     if (_y_ != 0 && _y_ != 0.5) _f_ |= FRPF_LOW;			\
     if (((r) >> _f_)&1) {						\
       for (;;) {							\
	 _w_ = (_w_ + 1)&MASK32;					\
	 if (_w_ || !_i_) break;					\
	 _z_->frac[_i_] = 0; _w_ = _z_->frac[--_i_];			\
       }								\
       if (!_w_) { _z_->exp++; _w_ = B31; }				\
       _z_->frac[_i_] = w;						\
     }									\
									\
     /* The result is not exact. */					\
     (rc) |= FLTERR_INEXACT;						\
   } while (0)
#endif

#ifdef FROB_NANS
#  define FROBNAN_DEC do {						\
     if (_z->f&FLTF_NANMASK) _z->f ^= FLTF_NANMASK;			\
   } while (0)
#else
#  define FROBNAN_DEC do ; while (0)
#endif

#define DECFLT(ty, TY, frexp, rc, z_out, x, r) do {			\
  struct floatbits *_z = (z_out);					\
  unsigned _rc = 0;							\
									\
  switch (TY##_FORMAT&(FLTFMT_ORGMASK | FLTFMT_TYPEMASK)) {		\
									\
    case FLTFMT_IEEE_F32: {						\
      unsigned _t[1];							\
      unsigned char *_x = (unsigned char *)&(x);			\
									\
      switch (TY##_FORMAT&FLTFMT_ENDMASK) {				\
	case FLTFMT_BE: _t[0] = LOAD32_B(_x); break;			\
	case FLTFMT_LE: _t[0] = LOAD32_L(_x); break;			\
	default: assert(!"unimplemented byte order"); break;		\
      }									\
      _rc |= fltfmt_decieee(&fltfmt_f32, _z, _t); FROBNAN_DEC;		\
    } break;								\
									\
    case FLTFMT_IEEE_F64: {						\
      unsigned _t[2];							\
      unsigned char *_x = (unsigned char *)&(x);			\
									\
      switch (TY##_FORMAT&FLTFMT_ENDMASK) {				\
	case FLTFMT_BE:							\
	  _t[0] = LOAD32_B(_x + 0); _t[1] = LOAD32_B(_x + 4);		\
	  break;							\
	case FLTFMT_LE:							\
	  _t[1] = LOAD32_L(_x + 0); _t[0] = LOAD32_L(_x + 4);		\
	  break;							\
	case FLTFMT_ARME:						\
	  _t[0] = LOAD32_L(_x + 0); _t[1] = LOAD32_L(_x + 4);		\
	  break;							\
	default: assert(!"unimplemented byte order"); break;		\
      }									\
      _rc |= fltfmt_decieee(&fltfmt_f64, _z, _t); FROBNAN_DEC;		\
    } break;								\
									\
    case FLTFMT_IEEE_F128: {						\
      unsigned _t[4];							\
      unsigned char *_x = (unsigned char *)&(x);			\
									\
      switch (TY##_FORMAT&FLTFMT_ENDMASK) {				\
	case FLTFMT_BE:							\
	  _t[0] = LOAD32_B(_x +  0); _t[1] = LOAD32_B(_x +  4);		\
	  _t[2] = LOAD32_B(_x +  8); _t[3] = LOAD32_B(_x + 12);		\
	  break;							\
	case FLTFMT_LE:							\
	  _t[3] = LOAD32_L(_x +  0); _t[2] = LOAD32_L(_x +  4);		\
	  _t[1] = LOAD32_L(_x +  8); _t[0] = LOAD32_L(_x + 12);		\
	  break;							\
	default: assert(!"unimplemented byte order"); break;		\
      }									\
      _rc |= fltfmt_decieee(&fltfmt_f128, _z, _t); FROBNAN_DEC;		\
    } break;								\
									\
    case FLTFMT_INTEL_F80: {						\
      unsigned _t[3];							\
      unsigned char *_x = (unsigned char *)&(x);			\
									\
      switch (TY##_FORMAT&FLTFMT_ENDMASK) {				\
	case FLTFMT_BE:							\
	  _t[0] = LOAD16_B(_x + 0);					\
	  _t[1] = LOAD32_B(_x + 2); _t[2] = LOAD32_B(_x + 6);		\
	  break;							\
	case FLTFMT_LE:							\
	  _t[2] = LOAD32_L(_x + 0); _t[1] = LOAD32_L(_x + 4);		\
	  _t[0] = LOAD16_L(_x + 8);					\
	  break;							\
	default: assert(!"unimplemented byte order"); break;		\
      }									\
      _rc |= fltfmt_decieee(&fltfmt_idblext80, _z, _t); FROBNAN_DEC;	\
    } break;								\
									\
    default: {								\
      ty _x = (x), _y;							\
      unsigned _i, _n, _f = 0;						\
      uint32 _t;							\
									\
      /* If the value looks negative then negate it and set the sign	\
       * flag.								\
       */								\
      if (NEGP(_x)) { _f |= FLTF_NEG; _x = -_x; }			\
									\
      /* Now for the case analysis.  Infinities and zero are		\
       * unproblematic.  NaNs can't be decoded exactly using the	\
       * portable machinery.						\
       */								\
      if (INFP(_x)) _f |= FLTF_INF;					\
      else if (_x == 0.0) _f |= FLTF_ZERO;				\
      else if (NANP(_x)) { _f |= FLTF_QNAN; _rc |= FLTERR_INEXACT; }	\
      else {								\
	/* A finite value.  Determine the number of fraction limbs	\
	 * we'll need based on the precision and radix and pull out	\
	 * 32-bit chunks one at a time.  This will be unproblematic	\
	 * for power-of-two radices, requiring at most shifting the	\
	 * significand left by a few bits, but inherently inexact (for	\
	 * the most part) for others.					\
	 */								\
									\
	_n = NFRAC(TY); fltfmt_allocfrac(_z, _n);			\
	_y = frexp(_x, &_z->exp);					\
	for (_i = 0; _i < _n; _i++)					\
	  { _y *= SH32; _t = _y; _y -= _t; _z->frac[_i] = _t; }		\
	CONVFIX(ty, _rc, _z, _y, _n, _f, (r));				\
      }									\
									\
      /* Done. */							\
      _z->f = _f;							\
    } break;								\
  }									\
									\
  /* Set the error code. */						\
  (rc) = _rc;								\
} while (0)

/* --- @fltfmt_decTY@ --- *
 *
 * Arguments:	@struct floatbits *z_out@ = storage for the result
 *		@ty x@ = value to decode
 *		@unsigned r@ = rounding mode
 *
 * Returns:	Error flags (@FLTERR_...@).
 *
 * Use:		Decode the native C floatingpoint value @x@ and store the
 *		result in @z_out@.
 *
 *		The @TY@ may be @flt@ to encode a @float@, @dbl@ to encode a
 *		@double@, or (on C99 implementations) @ldbl@ to encode a
 *		@long double@.
 *
 *		In detail, conversion is performed as follows.
 *
 *		  * If the implementation supports negative zeros and/or
 *		    infinity, then these are recognized and decoded.
 *
 *		  * If the input as a NaN, but the implementation cannot
 *		    usefully report NaN payloads, then the @FLTERR_INEXACT@
 *		    error bit is set and the decoded payload is left empty.
 *
 *		  * If the implementation's floating-point radix is not a
 *		    power of two, and @x@ is a nonzero finite value, then
 *		    @FLTERR_INEXACT@ error bit is set (unconditionally), and
 *		    the rounded value (according to the rounding mode @r@) is
 *		    stored in as many fraction words as necessary to identify
 *		    the original value uniquely.  If the radix is a power of
 *		    two, then the value is represented exactly.
 */

unsigned fltfmt_decflt(struct floatbits *z_out, float x, unsigned r)
{
  unsigned rc;

  DECFLT(double, FLT, frexp, rc, z_out, x, r);
  return (rc);
}

unsigned fltfmt_decdbl(struct floatbits *z_out, double x, unsigned r)
{
  unsigned rc;

  DECFLT(double, DBL, frexp, rc, z_out, x, r);
  return (rc);
}

#if __STDC_VERSION__ >= 199001
unsigned fltfmt_decldbl(struct floatbits *z_out, long double x, unsigned r)
{
  unsigned rc;

  DECFLT(long double, LDBL, frexpl, rc, z_out, x, r);
  return (rc);
}
#endif

/*----- That's all, folks -------------------------------------------------*/