Max 5 API Reference

jit.half.c

/*
 * Copyright 2001-2005 - Cycling '74
 * Derek Gerstmann - derek@cycling74.com
 *
 * Support for 16 bit floating point evaluation.  Based on the OpenEXR
 * C++ half datatype.  License information follows.
 *
 */

#include "jit.common.h"
#include "jit.half.h"


///////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
// Digital Ltd. LLC
// 
// All rights reserved.
// 
// 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.
// *       Neither the name of Industrial Light & Magic nor the names of
// its contributors may be used to endorse or promote products derived
// from this software without specific prior written permission. 
// 
// 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.
//
///////////////////////////////////////////////////////////////////////////

// --------------------------------------------------------------------------

// combined unsigned int & float union
typedef union _uif
{
    unsigned int i;
    float f;
} uif;

// --------------------------------------------------------------------------

static long jit_half_flut_ok = FALSE;   // init flag for half-to-float lut
static long jit_half_elut_ok = FALSE;   // init flag for float-to-half lut
static uif jit_half_flut[1 << 16];      // half-to-float lookup table
static half jit_half_elut[1 << 9];      // float-to-half exponent lookup table

// --------------------------------------------------------------------------

static void jit_half_elut_init();
static void jit_half_flut_init();
static half jit_bits_to_half(int i);
static unsigned int jit_half_to_bits(half y);
static float jit_half_overflow();

//---------------------------------------------------------------------------
//
// Implementation --
//
// Representation of a float:
//
//  We assume that a float, f, is an IEEE 754 single-precision
//  floating point number, whose bits are arranged as follows:
//
//      31 (msb)
//      | 
//      | 30     23
//      | |      | 
//      | |      | 22                    0 (lsb)
//      | |      | |                     |
//      X XXXXXXXX XXXXXXXXXXXXXXXXXXXXXXX
//
//      s e        m
//
//  S is the sign-bit, e is the exponent and m is the significand.
//
//  If e is between 1 and 254, f is a normalized number:
//
//              s    e-127
//      f = (-1)  * 2      * 1.m
//
//  If e is 0, and m is not zero, f is a denormalized number:
//
//              s    -126
//      f = (-1)  * 2      * 0.m
//
//  If e and m are both zero, f is zero:
//
//      f = 0.0
//
//  If e is 255, f is an "infinity" or "not a number" (NAN),
//  depending on whether m is zero or not.
//
//  Examples:
//
//      0 00000000 00000000000000000000000 = 0.0
//      0 01111110 00000000000000000000000 = 0.5
//      0 01111111 00000000000000000000000 = 1.0
//      0 10000000 00000000000000000000000 = 2.0
//      0 10000000 10000000000000000000000 = 3.0
//      1 10000101 11110000010000000000000 = -124.0625
//      0 11111111 00000000000000000000000 = +infinity
//      1 11111111 00000000000000000000000 = -infinity
//      0 11111111 10000000000000000000000 = NAN
//      1 11111111 11111111111111111111111 = NAN
//
// Representation of a half:
//
//  Here is the bit-layout for a half number, h:
//
//      15 (msb)
//      | 
//      | 14  10
//      | |   |
//      | |   | 9        0 (lsb)
//      | |   | |        |
//      X XXXXX XXXXXXXXXX
//
//      s e     m
//
//  S is the sign-bit, e is the exponent and m is the significand.
//
//  If e is between 1 and 30, h is a normalized number:
//
//              s    e-15
//      h = (-1)  * 2     * 1.m
//
//  If e is 0, and m is not zero, h is a denormalized number:
//
//              S    -14
//      h = (-1)  * 2     * 0.m
//
//  If e and m are both zero, h is zero:
//
//      h = 0.0
//
//  If e is 31, h is an "infinity" or "not a number" (NAN),
//  depending on whether m is zero or not.
//
//  Examples:
//
//      0 00000 0000000000 = 0.0
//      0 01110 0000000000 = 0.5
//      0 01111 0000000000 = 1.0
//      0 10000 0000000000 = 2.0
//      0 10000 1000000000 = 3.0
//      1 10101 1111000001 = -124.0625
//      0 11111 0000000000 = +infinity
//      1 11111 0000000000 = -infinity
//      0 11111 1000000000 = NAN
//      1 11111 1111111111 = NAN
//
// Conversion:
//
//  Converting from a float to a half requires some non-trivial bit
//  manipulations.  In some cases, this makes conversion relatively
//  slow, but the most common case is accelerated via table lookups.
//
//  Converting back from a half to a float is easier because we don't
//  have to do any rounding.  In addition, there are only 65536
//  different half numbers; we can convert each of those numbers once
//  and store the results in a table.  Later, all conversions can be
//  done using only simple table lookups.
//
//---------------------------------------------------------------------------

// --------------------------------------------------------------------------
// Convert a float to half
// --------------------------------------------------------------------------
half jit_float_to_half(float f)
{
    uif x;
    half h;
    register int e;
    if (f == 0)
    {
        // common special case - zero.
        // for speed, we don't preserve the zero's sign.
        h = 0;
    }
    else
    {
        // We extract the combined sign and exponent, e, from our
        // floating-point number, f.  Then we convert e to the sign
        // and exponent of the half number via a table lookup.
        //
        // For the most common case, where a normalized half is produced,
        // the table lookup returns a non-zero value; in this case, all
        // we have to do, is round f's significand to 10 bits and combine
        // the result with e.
        //
        // For all other cases (overflow, zeroes, denormalized numbers
        // resulting from underflow, infinities and NANs), the table
        // lookup returns zero, and we call a longer, non-inline function
        // to do the float-to-half conversion.
        x.f = f;
        e = (x.i >> 23) & 0x000001ff;
        e = jit_half_elut[e];

        if (e)
        {
            // Simple case - round the significand and
            // combine it with the sign and exponent.
            h = e + (((x.i & 0x007fffff) + 0x00001000) >> 13);
        }
        else
        {
            // Difficult case - call a function.
            h = jit_bits_to_half(x.i);
        }
    }
    return h;
}

// --------------------------------------------------------------------------
// Round to n-bit precision
// --------------------------------------------------------------------------
half jit_half_round(half h, unsigned int n)
{

    unsigned short s, e;

    // Parameter check.
    if (n >= 10)
        return h;

    // Disassemble h into the sign, s,
    // and the combined exponent and significand, e.
    s = h & 0x8000;
    e = h & 0x7fff;

    // Round the exponent and significand to the nearest value
    // where ones occur only in the (10-n) most significant bits.
    // Note that the exponent adjusts automatically if rounding
    // up causes the significand to overflow.
    e >>= 9 - n;
    e  += e & 1;
    e <<= 9 - n;

    // Check for exponent overflow.
    if (e >= 0x7c00)
    {
        // Overflow occurred -- truncate instead of rounding.
        e = h;
        e >>= 10 - n;
        e <<= 10 - n;
    }

    // Put the original sign bit back.
    h = s | e;
    return h;
}

// --------------------------------------------------------------------------
// Convert a half-to-float
// --------------------------------------------------------------------------
float jit_half_to_float(half y)
{
    return jit_half_flut[y].f;
/*
    // use the following without a lookup table
    unsigned int r = jit_half_to_bits(y);
    return *(float*)(&r);
*/
}

// --------------------------------------------------------------------------
// Arithmetic operations
// --------------------------------------------------------------------------
half jit_half_negate(half h)
{
    return h ^ 0x8000;
}

half jit_half_add(half h, half a)
{
    return jit_half_to_float(jit_half_to_float(h) + jit_half_to_float (a));
}

half jit_half_add_float(half h, float a)
{
    return jit_half_to_float(jit_half_to_float(h) + a);
}

half jit_half_sub(half h, half a)
{
    return jit_float_to_half(jit_half_to_float(h) - jit_half_to_float(a));
}

half jit_half_sub_float(half h, float a)
{
    return jit_float_to_half(jit_half_to_float(h) - a);
}

half jit_half_mul(half h, half a)
{
    return jit_float_to_half(jit_half_to_float(h) * jit_half_to_float(a));
}

half jit_half_mul_float(half h, float a)
{
    return jit_float_to_half(jit_half_to_float(h) * a);
}

half jit_half_div(half h, half a)
{
    return jit_float_to_half(jit_half_to_float(h) / jit_half_to_float(a));
}

half jit_half_div_float(half h, float a)
{
    return jit_float_to_half(jit_half_to_float(h) / a);
}

long jit_half_is_finite(half h)
{
    unsigned short e = (h >> 10) & 0x001f;
    return e < 31;
}

long jit_half_is_normalized(half h)
{
    unsigned short e = (h >> 10) & 0x001f;
    return e > 0 && e < 31;
}

long jit_half_is_denormalized(half h)
{
    unsigned short e = (h >> 10) & 0x001f;
    unsigned short m =  h & 0x3ff;
    return e == 0 && m != 0;
}

long jit_half_is_zero(half h)
{
    return (h & 0x7fff) == 0;
}

long jit_half_is_negative(half h)
{
    return (h & 0x8000) != 0;
}


long jit_half_is_nan(half h)
{
    unsigned short e = (h >> 10) & 0x001f;
    unsigned short m =  h & 0x3ff;
    return e == 31 && m != 0;
}

long jit_half_is_inf(half h)
{
    unsigned short e = (h >> 10) & 0x001f;
    unsigned short m =  h & 0x3ff;
    return e == 31 && m == 0;
}

// --------------------------------------------------------------------------
// Initialize the lookup tables for half/float conversion
// --------------------------------------------------------------------------
void jit_half_init()
{
    if(!jit_half_elut_ok)
        jit_half_elut_init();
    jit_half_elut_ok = TRUE;

    if(!jit_half_flut_ok)
        jit_half_flut_init();
    jit_half_flut_ok = TRUE;
}

// --------------------------------------------------------------------------
// Compute a lookup table for float-to-half conversion.
//
// When indexed with the combined sign and exponent of
// a float, the table either returns the combined sign
// and exponent of the corresponding half, or zero if
// the corresponding half may not be normalized (zero,
// denormalized, overflow).
// --------------------------------------------------------------------------
static void jit_half_elut_init()
{
    int i, e;
    for (i = 0; i < 0x100; i++)
    {
        e = (i & 0x0ff) - (127 - 15);

        if (e <= 0 || e >= 30)
        {
            // special case
            jit_half_elut[i]         = 0;
            jit_half_elut[i | 0x100] = 0;
        }
        else
        {
            // common case - normalized half, no exponent overflow possible
            jit_half_elut[i]         =  (e << 10);
            jit_half_elut[i | 0x100] = ((e << 10) | 0x8000);
        }
    }
}

// --------------------------------------------------------------------------
// Compute a lookup table for half-to-float conversion.
// --------------------------------------------------------------------------
static void jit_half_flut_init()
{
    int i;
    uif x;
    unsigned int ui;
    for(i = 0; i < (1 << 16); i++)
    {
        ui = jit_half_to_bits(i);
        x.i = ui;
        jit_half_flut[i] = x;
    }
}

// --------------------------------------------------------------------------
// Compose a half from a bitwise representation
// --------------------------------------------------------------------------
static half jit_bits_to_half(int i)
{
    // Our floating point number, f, is represented by the bit
    // pattern in integer i.  Disassemble that bit pattern into
    // the sign, s, the exponent, e, and the significand, m.
    // Shift s into the position where it will go in in the
    // resulting half number.
    // Adjust e, accounting for the different exponent bias
    // of float and half (127 versus 15).

    register int s =  (i >> 16) & 0x00008000;
    register int e = ((i >> 23) & 0x000000ff) - (127 - 15);
    register int m =   i        & 0x007fffff;

    // Now reassemble s, e and m into a half:
    if (e <= 0)
    {
        if (e < -10)
        {
            // E is less than -10.  The absolute value of f is
            // less than HALF_MIN (f may be a small normalized
            // float, a denormalized float or a zero).
            //
            // We convert f to a half zero.
            return 0;
        }

        // E is between -10 and 0.  F is a normalized float,
        // whose magnitude is less than HALF_NRM_MIN.
        //
        // We convert f to a denormalized half.
        m = (m | 0x00800000) >> (1 - e);

        // Round to nearest, round "0.5" up.
        //
        // Rounding may cause the significand to overflow and make
        // our number normalized.  Because of the way a half's bits
        // are laid out, we don't have to treat this case separately;
        // the code below will handle it correctly.
        if (m &  0x00001000)
            m += 0x00002000;

        // Assemble the half from s, e (zero) and m.
        return s | (m >> 13);
    }
    else if (e == 0xff - (127 - 15))
    {
        if (m == 0)
        {
            // F is an infinity; convert f to a half
            // infinity with the same sign as f.
            return s | 0x7c00;
        }
        else
        {
            // F is a NAN; we produce a half NAN that preserves
            // the sign bit and the 10 leftmost bits of the
            // significand of f, with one exception: If the 10
            // leftmost bits are all zero, the NAN would turn 
            // into an infinity, so we have to set at least one
            // bit in the significand.
            m >>= 13;
            return s | 0x7c00 | m | (m == 0);
        }
    }
    else
    {
        // E is greater than zero.  F is a normalized float.
        // We try to convert f to a normalized half.

        // Round to nearest, round "0.5" up
        if (m &  0x00001000)
        {
            m += 0x00002000;

            if (m & 0x00800000)
            {
            m =  0;     // overflow in significand,
            e += 1;     // adjust exponent
            }
        }

        // Handle exponent overflow
        if (e > 30)
        {   
            // cause a hardware floating point overflow
            jit_half_overflow();


            // if this returns, the half becomes an
            // infinity with the same sign as f.
            return s | 0x7c00;              

        }               

        // Assemble the half from s, e and m.
        return s | (e << 10) | (m >> 13);
    }
}

// --------------------------------------------------------------------------
// Decompose a half to a bitwise representation
// --------------------------------------------------------------------------
static unsigned int jit_half_to_bits(half y)
{
    register int s = (y >> 15) & 0x00000001;
    register int e = (y >> 10) & 0x0000001f;
    register int m =  y        & 0x000003ff;

    if (e == 0)
    {
        if (m == 0) // Plus or minus zero
        {
            return s << 31;
        }
        else // Denormalized number -- renormalize it
        {
            while (!(m & 0x00000400))
            {
                m <<= 1;
                e -=  1;
            }

            e += 1;
            m &= ~0x00000400;
        }
    }
    else if (e == 31)
    {
        if (m == 0) // Inf
        {
            return (s << 31) | 0x7f800000;
        }
        else // NaN
        {
            return (s << 31) | 0x7f800000 | (m << 13);
        }
    }

    e = e + (127 - 15);
    m = m << 13;

    return (s << 31) | (e << 23) | m;
}

// --------------------------------------------------------------------------
// Overflow handler for float-to-half conversion which generates a 
// hardware floating-point overflow, which may be trapped by the OS.
// --------------------------------------------------------------------------
static float jit_half_overflow()
{
    int i;
    volatile float f = 1e10;

    for (i = 0; i < 10; i++)    
        f *= f;             // this will overflow before
                            // the for�loop terminates
    return f;
}

// --------------------------------------------------------------------------

Copyright © 2008, Cycling '74