Add utilities to convert between half and single precisionhalf-precision

Derived from BSD licensed code written by James Tursa, downloaded from
here:

   http://www.mathworks.com/matlabcentral/fileexchange/23173

2 files changed, 307 insertions, 0 deletions

diff --git a/pixman/pixman-private.h b/pixman/pixman-private.h
index 664260b9..a650918a 100644
--- a/pixman/pixman-private.h
+++ b/pixman/pixman-private.h
@@ -572,6 +572,12 @@ _pixman_choose_implementation (void);
 uint32_t *
 _pixman_iter_get_scanline_noop (pixman_iter_t *iter, const uint32_t *mask);
 
+void
+_pixman_float_to_half (uint16_t *target, float *source, int n_elements);
+
+void
+_pixman_half_to_float (float *target, uint16_t *source, int n_elements);
+
 /* These "formats" all have depth 0, so they
  * will never clash with any real ones
  */
diff --git a/pixman/pixman-utils.c b/pixman/pixman-utils.c
index cb4e6219..ca6bd9a7 100644
--- a/pixman/pixman-utils.c
+++ b/pixman/pixman-utils.c
@@ -242,6 +242,307 @@ pixman_region32_copy_from_region16 (pixman_region32_t *dst,
     return retval;
 }
 
+/* Half <-> Float conversion */
+
+/******************************************************************************
+ *
+ * Filename:    ieeehalfprecision.c
+ * Programmer:  James Tursa
+ * Version:     1.0
+ * Date:        March 3, 2009
+ * Copyright:   (c) 2009 by James Tursa, All Rights Reserved
+ *
+ *  This code uses the BSD License:
+ *
+ *  Redistribution and use in source and binary forms, with or without
+ *  modification, are permitted provided that the following conditions are
+ *  met:
+ *
+ *     * Redistributions of source code must retain the above copyright
+ *       notice, this list of conditions and the following disclaimer.
+ *     * Redistributions in binary form must reproduce the above copyright
+ *       notice, this list of conditions and the following disclaimer in
+ *       the documentation and/or other materials provided with the distribution
+ *
+ *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ *  POSSIBILITY OF SUCH DAMAGE.
+ *
+ * This file contains C code to convert between IEEE double, single, and half
+ * precision floating point formats. The intended use is for standalone C code
+ * that does not rely on MATLAB mex.h. The bit pattern for the half precision
+ * floating point format is stored in a 16-bit unsigned int variable. The half
+ * precision bit pattern definition is:
+ *
+ * 1 bit sign bit
+ * 5 bits exponent, biased by 15
+ * 10 bits mantissa, hidden leading bit, normalized to 1.0
+ *
+ * Special floating point bit patterns recognized and supported:
+ *
+ * All exponent bits zero:
+ * - If all mantissa bits are zero, then number is zero (possibly signed)
+ * - Otherwise, number is a denormalized bit pattern
+ *
+ * All exponent bits set to 1:
+ * - If all mantissa bits are zero, then number is +Infinity or -Infinity
+ * - Otherwise, number is NaN (Not a Number)
+ *
+ * For the denormalized cases, note that 2^(-24) is the smallest number that
+ * can be represented in half precision exactly. 2^(-25) will convert to
+ * 2^(-24) because of the rounding algorithm used, and 2^(-26) is too small
+ * and underflows to zero.
+ *
+ *****************************************************************************/
+
+/*
+ * Modifications to this code by Soren Sandmann:
+ *  - delete checkieee bits
+ *  - delete double routines
+ *  - don't allow NULL source/target
+ *  - make routines return void instead of int
+ *  - convert C++ style comments to C comments
+ *  - reformatting
+ *
+ * FIXME: Consider preserving as much of the NaN payload as possible.
+ * Currently, the mantissa is just set to 0 when a NaN is converted.
+ * This is what Ivy Bridge is supposed to do - see the Intel AVX manual
+ * for details.
+ */
+
+#include <string.h>
+#include <stdint.h>
+
+/*-----------------------------------------------------------------------------
+ *
+ * Routine:  single_to_half
+ *
+ * Input:    source = Address of 32-bit floating point data to convert
+ *
+ * Output:   target = Address of 16-bit data to hold output
+ *
+ * Programmer:  James Tursa
+ *
+ *-----------------------------------------------------------------------------
+ */
+void
+pixman_float_to_half (uint16_t *target, float *source, int n_elements)
+{
+    uint16_t * hp = (uint16_t *) target;
+    uint32_t * xp = (uint32_t *) source;
+    uint16_t   hs, he, hm;
+    uint32_t   x, xs, xe, xm;
+    int        hes;
+
+    while (n_elements--)
+    {
+	x = *xp++;
+	if ((x & 0x7FFFFFFFu) == 0)
+	{
+	    /* Signed zero */
+	    *hp++ = (uint16_t) (x >> 16);
+	}
+	else
+	{
+	    xs = x & 0x80000000u;  /* Pick off sign bit */
+	    xe = x & 0x7F800000u;  /* Pick off exponent bits */
+	    xm = x & 0x007FFFFFu;  /* Pick off mantissa bits */
+
+	    if (xe == 0)
+	    {
+		/* Denormal will underflow, return a signed zero */
+
+		*hp++ = (uint16_t) (xs >> 16);
+	    }
+	    else if (xe == 0x7F800000u) /* Inf or NaN (all exponent bits set */
+	    {
+		if (xm == 0)
+		{
+		    /* Signed Inf */
+		    *hp++ = (uint16_t)((xs >> 16) | 0x7C00u);
+		}
+		else
+		{
+		    /* NaN, only 1st mantissa bit set */
+		    *hp++ = (uint16_t)0xFE00u;
+		}
+	    }
+	    else /* Normalized number */
+	    {
+		/* Sign bit */
+		hs = (uint16_t) (xs >> 16);
+
+		/* Exponent unbias the single, then bias the halfp */
+		hes = ((int)(xe >> 23)) - 127 + 15;
+
+		if (hes >= 0x1F)
+		{
+		    /* Overflow => signed Inf */
+		    *hp++ = (uint16_t) ((xs >> 16) | 0x7C00u);
+		}
+		else if (hes <= 0)
+		{
+		    /* Underflow */
+		    if ((14 - hes) > 24)
+		    {
+			/* Mantissa shifted all the way off
+			 * no rounding possibility
+			 */
+			hm = (uint16_t) 0u;  /* Set mantissa to zero */
+		    }
+		    else
+		    {
+			/* Add the hidden leading bit */
+			xm |= 0x00800000u;
+
+			hm = (uint16_t) (xm >> (14 - hes)); /* Mantissa */
+
+			/* Check for rounding */
+			if ((xm >> (13 - hes)) & 0x00000001u)
+			{
+			    /* Round, might overflow into exp bit,
+			     * but this is OK
+			     */
+			    hm += (uint16_t) 1u;
+			}
+		    }
+
+		    /* Combine sign and mantissa, biased exponent is zero */
+		    *hp++ = (hs | hm);
+		}
+		else
+		{
+		    /* Exponent */
+		    he = (uint16_t) (hes << 10);
+
+		    /* Mantissa */
+		    hm = (uint16_t) (xm >> 13);
+
+		    /* Check for rounding */
+		    if (xm & 0x00001000u)
+		    {
+			/* Round, might overflow to inf, but this is OK */
+			*hp++ = (hs | he | hm) + (uint16_t) 1u;
+		    }
+		    else
+		    {
+			/* No rounding */
+			*hp++ = (hs | he | hm);
+		    }
+		}
+	    }
+	}
+    }
+}
+
+/*-----------------------------------------------------------------------------
+ *
+ * Routine:  half_to_singles
+ *
+ * Input:  source = address of 16-bit data to convert
+ *         n_elements  = Number of values at that address to convert
+ *
+ * Output: target = Address of 32-bit floating point data to hold output
+ *
+ * Programmer:  James Tursa
+ *
+ *-----------------------------------------------------------------------------
+ */
+void
+pixman_half_to_float (float *target, uint16_t *source, int n_elements)
+{
+    uint16_t *hp = (uint16_t *) source;
+    uint32_t *xp = (uint32_t *) target;
+    uint16_t  h, hs, he, hm;
+    uint32_t  xs, xe, xm;
+    int32_t   xes;
+    int       e;
+
+    while (n_elements--)
+    {
+	h = *hp++;
+
+	if ((h & 0x7FFFu) == 0)    /* Signed zero */
+	{
+	    *xp++ = ((uint32_t) h) << 16;  /* Return the signed zero */
+	}
+	else     /* Not zero */
+	{
+	    hs = h & 0x8000u;  /* Pick off sign bit */
+	    he = h & 0x7C00u;  /* Pick off exponent bits */
+	    hm = h & 0x03FFu;  /* Pick off mantissa bits */
+
+	    if (he == 0) /* Denormal will convert to normalized */
+	    {
+		/* Figure out how much extra to adjust the exponent */
+		e = -1;
+		do
+		{
+		    e++;
+		    hm <<= 1;
+
+		    /* Shift until leading bit overflows into exponent bit */
+		}
+		while ((hm & 0x0400u) == 0);
+
+		/* Sign bit */
+		xs = ((uint32_t) hs) << 16; /* Sign bit */
+
+		/* Exponent unbias the halfp, then bias the single */
+		xes = ((int32_t) (he >> 10)) - 15 + 127 - e;
+
+		/* Exponent */
+		xe = (uint32_t) (xes << 23);
+
+		/* Mantissa */
+		xm = ((uint32_t) (hm & 0x03FFu)) << 13;
+
+		/* Combine sign, exponent and mantissa */
+		*xp++ = (xs | xe | xm); 
+	    }
+	    else if (he == 0x7C00u)  /* Inf or NaN (all exponent bits set) */
+	    {
+		/* If mantissa is zero ... */
+		if (hm == 0)
+		{
+		    /* Signed Inf */
+		    *xp++ = (((uint32_t) hs) << 16) | ((uint32_t) 0x7F800000u);
+		}
+		else
+		{
+		    /* NaN, only 1st mantissa bit set */
+		    *xp++ = (uint32_t) 0xFFC00000u;
+		}
+	    }
+	    else /* Normalized number */
+	    {
+		/* Sign bit */
+		xs = ((uint32_t) hs) << 16;
+
+		/* Unbias the halfp, then bias the single */
+		xes = ((int32_t) (he >> 10)) - 15 + 127;
+
+		/* Exponent */
+		xe = (uint32_t) (xes << 23);             /* Exponent */
+
+		/* Mantissa */
+		xm = ((uint32_t) hm) << 13;              /* Mantissa */
+
+		/* Combine sign, exponent and mantissa */
+		*xp++ = (xs | xe | xm);
+	    }
+	}
+    }
+}
+
 #ifdef DEBUG
 
 void