pack (const unsigned short s[16],
      unsigned char b[14],
      bool optFlatFields,
      bool exactMax)
{
    //
    // Pack a block of 4 by 4 16-bit pixels (32 bytes) into
    // either 14 or 3 bytes.
    //

    //
    // Integers s[0] ... s[15] represent floating-point numbers
    // in what is essentially a sign-magnitude format.  Convert
    // s[0] .. s[15] into a new set of integers, t[0] ... t[15],
    // such that if t[i] is greater than t[j], the floating-point
    // number that corresponds to s[i] is always greater than
    // the floating-point number that corresponds to s[j].
    //
    // Also, replace any bit patterns that represent NaNs or
    // infinities with bit patterns that represent floating-point
    // zeroes.
    //
    //	bit pattern	floating-point		bit pattern
    //	in s[i]		value			in t[i]
    //
    //  0x7fff		NAN			0x8000
    //  0x7ffe		NAN			0x8000
    //	  ...					  ...
    //  0x7c01		NAN			0x8000
    //  0x7c00		+infinity		0x8000
    //  0x7bff		+HALF_MAX		0xfbff
    //  0x7bfe					0xfbfe
    //  0x7bfd					0xfbfd
    //	  ...					  ...
    //  0x0002		+2 * HALF_MIN		0x8002
    //  0x0001		+HALF_MIN		0x8001
    //  0x0000		+0.0			0x8000
    //  0x8000		-0.0			0x7fff
    //  0x8001		-HALF_MIN		0x7ffe
    //  0x8002		-2 * HALF_MIN		0x7ffd
    //	  ...					  ...
    //  0xfbfd					0x0f02
    //  0xfbfe					0x0401
    //  0xfbff		-HALF_MAX		0x0400
    //  0xfc00		-infinity		0x8000
    //  0xfc01		NAN			0x8000
    //	  ...					  ...
    //  0xfffe		NAN			0x8000
    //  0xffff		NAN			0x8000
    //

    unsigned short t[16];

    for (int i = 0; i < 16; ++i)
    {
	if ((s[i] & 0x7c00) == 0x7c00)
	    t[i] = 0x8000;
	else if (s[i] & 0x8000)
	    t[i] = ~s[i];
	else
	    t[i] = s[i] | 0x8000;
    }
    
    //
    // Find the maximum, tMax, of t[0] ... t[15].
    //

    unsigned short tMax = 0;

    for (int i = 0; i < 16; ++i)
	if (tMax < t[i])
	    tMax = t[i];

    //
    // Compute a set of running differences, r[0] ... r[14]:
    // Find a shift value such that after rounding off the
    // rightmost bits and shifting all differenes are between
    // -32 and +31.  Then bias the differences so that they
    // end up between 0 and 63.
    //

    int shift = -1;
    int d[16];
    int r[15];
    int rMin;
    int rMax;

    const int bias = 0x20;

    do
    {
        shift += 1;

        //
        // Compute absolute differences, d[0] ... d[15],
        // between tMax and t[0] ... t[15].
        //
        // Shift and round the absolute differences.
        //

        for (int i = 0; i < 16; ++i)
            d[i] = shiftAndRound (tMax - t[i], shift);

        //
        // Convert d[0] .. d[15] into running differences
        //

        r[ 0] = d[ 0] - d[ 4] + bias;
        r[ 1] = d[ 4] - d[ 8] + bias;
        r[ 2] = d[ 8] - d[12] + bias;

        r[ 3] = d[ 0] - d[ 1] + bias;
        r[ 4] = d[ 4] - d[ 5] + bias;
        r[ 5] = d[ 8] - d[ 9] + bias;
        r[ 6] = d[12] - d[13] + bias;

        r[ 7] = d[ 1] - d[ 2] + bias;
        r[ 8] = d[ 5] - d[ 6] + bias;
        r[ 9] = d[ 9] - d[10] + bias;
        r[10] = d[13] - d[14] + bias;

        r[11] = d[ 2] - d[ 3] + bias;
        r[12] = d[ 6] - d[ 7] + bias;
        r[13] = d[10] - d[11] + bias;
        r[14] = d[14] - d[15] + bias;

        rMin = r[0];
        rMax = r[0];

        for (int i = 1; i < 15; ++i)
        {
            if (rMin > r[i])
                rMin = r[i];

            if (rMax < r[i])
                rMax = r[i];
        }
    }
    while (rMin < 0 || rMax > 0x3f);

    if (rMin == bias && rMax == bias && optFlatFields)
    {
        //
        // Special case - all pixels have the same value.
        // We encode this in 3 instead of 14 bytes by
        // storing the value 0xfc in the third output byte,
        // which cannot occur in the 14-byte encoding.
        //

        b[0] = (t[0] >> 8);
        b[1] = (unsigned char) t[0];
        b[2] = 0xfc;

        return 3;
    }

    if (exactMax)
    {
        //
        // Adjust t[0] so that the pixel whose value is equal
        // to tMax gets represented as accurately as possible.
        //

        t[0] = tMax - (d[0] << shift);
    }

    //
    // Pack t[0], shift and r[0] ... r[14] into 14 bytes:
    //

    b[ 0] = (t[0] >> 8);
    b[ 1] = (unsigned char) t[0];

    b[ 2] = (unsigned char) ((shift << 2) | (r[ 0] >> 4));
    b[ 3] = (unsigned char) ((r[ 0] << 4) | (r[ 1] >> 2));
    b[ 4] = (unsigned char) ((r[ 1] << 6) |  r[ 2]      );

    b[ 5] = (unsigned char) ((r[ 3] << 2) | (r[ 4] >> 4));
    b[ 6] = (unsigned char) ((r[ 4] << 4) | (r[ 5] >> 2));
    b[ 7] = (unsigned char) ((r[ 5] << 6) |  r[ 6]      );

    b[ 8] = (unsigned char) ((r[ 7] << 2) | (r[ 8] >> 4));
    b[ 9] = (unsigned char) ((r[ 8] << 4) | (r[ 9] >> 2));
    b[10] = (unsigned char) ((r[ 9] << 6) |  r[10]      );

    b[11] = (unsigned char) ((r[11] << 2) | (r[12] >> 4));
    b[12] = (unsigned char) ((r[12] << 4) | (r[13] >> 2));
    b[13] = (unsigned char) ((r[13] << 6) |  r[14]      );

    return 14;
}