#
# how many numbers can a minifloat represent between 0 and 1?
# let's find out.
#
# a minifloat is a 16-bit s10e5 float:
# S Exp.  Mantissa
# x xxxxx xxxxxxxxxx
#
# Exponent is biased by -15.  e.g. exponent of 0 is -15, exponent of 14 is -1.
# If Exponent is 0, mantissa is denormalized, no automatic leading 1.
# Otherwise it has an automatic leading 1.
#
# 0 00000 0000000000 zero
#
# http://home.earthlink.net/~mrob/pub/math/floatformats.html#minifloat

"""
http://oss.sgi.com/projects/ogl-sample/registry/ARB/half_float_pixel.txt
 A 16-bit floating-point number has a 1-bit sign (S), a 5-bit
    exponent (E), and a 10-bit mantissa (M).  The value of a 16-bit
    floating-point number is determined by the following:

        (-1)^S * 0.0,                        if E == 0 and M == 0,
        (-1)^S * 2^-14 * (M / 2^10),         if E == 0 and M != 0,
        (-1)^S * 2^(E-15) * (1 + M/2^10),    if 0 < E < 31,
        (-1)^S * INF,                        if E == 31 and M == 0, or
        NaN,                                 if E == 31 and M != 0,

    where

        S = floor((N mod 65536) / 32768),
        E = floor((N mod 32768) / 1024), and
        M = N mod 1024.

"""

import sys


def split(i):
    sign = (i % 65536) / 32768
    exponent = (i % 32768) / 1024
    mantissa = i % 1024
    return (sign, exponent, mantissa)

twoToTheNegativeFourteenth = pow(2, -14)

def floatify(i):
    global twoToTheNegativeFourteenth
    global treatDenormalizedAsZero
    (sign, exponent, mantissa) = split(i)
    f = 0
    sign = 1 if (sign == 0) else -1

    if exponent == 0:
        f = sign * twoToTheNegativeFourteenth * (float(mantissa) / (2**10))
    elif exponent == 31:
        return sign * 5555555555555555555555555555555555555555555555.0 # too lazy to remember how to generate a real "infinity"
    else:
        f = sign * pow(2, exponent - 15) * (1 + (float(mantissa) / (2**10)))
    # print "f " + str(f)
    return f

if 0:
    floatify(1)
    floatify(1023)
    floatify(1024)
    print "float of 32767 is " + str(floatify(32767)) # infinity
    print "float of " + str(32767 & ~1024) + " is " + str(floatify(32767 & ~1024)) # largest number


numbersWithDenormalized = []
numbersWithoutDenormalized = []

for i in range(0, 65536):
    (sign, exponent, mantissa) = split(i)
    f = floatify(i)
    numbersWithDenormalized.append(f)
    if exponent == 0:
        numbersWithoutDenormalized.append(0.0)
    else:
        numbersWithoutDenormalized.append(f)


numbersWithDenormalized.sort()
numbersWithoutDenormalized.sort()

zeroToOneWithDenormalized = []
zeroToOneWithoutDenormalized = []

last = -11111 # sentinel value
for i in numbersWithDenormalized:
    if i == last:
        continue
    last = i
    if (i < 0) or (i > 1):
        continue
    zeroToOneWithDenormalized.append(i)

last = -11111 # sentinel value
for i in numbersWithoutDenormalized:
    if i == last:
        continue
    last = i
    if (i < 0) or (i > 1):
        continue
    zeroToOneWithoutDenormalized.append(i)

print "       Including all denormalized numbers, there are " + str(len(zeroToOneWithDenormalized)) + " m10e5 minifloats in the range (0, 1)."
print "Rounding all denormalized numbers to zero, there are " + str(len(zeroToOneWithoutDenormalized)) + " m10e5 minifloats in the range (0, 1)."

assert(zeroToOneWithDenormalized[0] == 0)
assert(zeroToOneWithoutDenormalized[0] == 0)

# if you'd like to see all the numbers, with the denormalized numbers starred,
# change this to "if 1:"
if 0:
    for i in zeroToOneWithDenormalized:
        starred = " *" if ((i > 0) and (i < zeroToOneWithoutDenormalized[1])) else ""
        print str(i) + starred

for exponent in range(11, 13):
    divisor = 2 ** exponent
    floatDivisor = float(divisor)
    failCount = 0
    for i in range(0, divisor + 1):
        f = i / floatDivisor
        if not f in zeroToOneWithoutDenormalized:
            failCount += 1
    if failCount == 0:
        print "zeroToOneWithoutDenormalized contains all values from 0 to 1 stepping by " + str(divisor) + " (precision of " + str(exponent) + " bits)."
    else:
        print " zeroToOneWithoutDenormalized lacked " + str(failCount) + " values from 0 to 1 stepping by " + str(divisor) + " (precision of " + str(exponent) + " bits)."