Here is code that calculates the number of values representable in float in all finite ranges. It expects IEEE-754 arithmetic. I adapted it from my previous C++ answer.
This has two implementations of converting a floating-point number to its encoding (one by copying the bits, one by mathematically manipulating it). After that, the distance calculation is fairly simple (negative values have to be adjusted, and then the distance is simply a subtraction).
#include <float.h>
#include <inttypes.h>
#include <limits.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <tgmath.h>
/* Define a value with only the high bit of a uint32_t set. This is also the
encoding of floating-point -0.
*/
static const uint32_t HighBit = UINT32_MAX ^ UINT32_MAX>>1;
// Return the encoding of a floating-point number by copying its bits.
static uint32_t EncodingBits(float x)
{
uint32_t result;
memcpy(&result, &x, sizeof result);
return result;
}
// Return the encoding of a floating-point number by using math.
static uint32_t EncodingMath(float x)
{
static const int SignificandBits = FLT_MANT_DIG;
static const int MinimumExponent = FLT_MIN_EXP;
// Encode the high bit.
uint32_t result = signbit(x) ? HighBit : 0;
// If the value is zero, the remaining bits are zero, so we are done.
if (x == 0) return result;
/* The C library provides a little-known routine to split a floating-point
number into a significand and an exponent. Note that this produces a
normalized significand, not the actual significand encoding. Notably,
it brings significands of subnormals up to at least 1/2. We will
adjust for that below. Also, this routine normalizes to [1/2, 1),
whereas IEEE 754 is usually expressed with [1, 2), but that does not
bother us.
*/
int xe;
float xf = frexp(fabs(x), &xe);
// Test whether the number is subnormal.
if (xe < MinimumExponent)
{
/* For a subnormal value, the exponent encoding is zero, so we only
have to insert the significand bits. This scales the significand
so that its low bit is scaled to the 1 position and then inserts it
into the encoding.
*/
result |= (uint32_t) ldexp(xf, xe - MinimumExponent + SignificandBits);
}
else
{
/* For a normal value, the significand is encoded without its leading
bit. So we subtract .5 to remove that bit and then scale the
significand so its low bit is scaled to the 1 position.
*/
result |= (uint32_t) ldexp(xf - .5, SignificandBits);
/* The exponent is encoded with a bias of (in C++'s terminology)
MinimumExponent - 1. So we subtract that to get the exponent
encoding and then shift it to the position of the exponent field.
Then we insert it into the encoding.
*/
result |= ((uint32_t) xe - MinimumExponent + 1) << (SignificandBits-1);
}
return result;
}
/* Return the encoding of a floating-point number. For illustration, we
get the encoding with two different methods and compare the results.
*/
static uint32_t Encoding(float x)
{
uint32_t xb = EncodingBits(x);
uint32_t xm = EncodingMath(x);
if (xb != xm)
{
fprintf(stderr, "Internal error encoding %.99g.\n", x);
fprintf(stderr, "\tEncodingBits says %#" PRIx32 ".\n", xb);
fprintf(stderr, "\tEncodingMath says %#" PRIx32 ".\n", xm);
exit(EXIT_FAILURE);
}
return xb;
}
/* Return the distance from a to b as the number of values representable in
float from one to the other. b must be greater than or equal to a. 0 is
counted only once.
*/
static uint32_t Distance(float a, float b)
{
uint32_t ae = Encoding(a);
uint32_t be = Encoding(b);
/* For represented values from +0 to infinity, the IEEE 754 binary
floating-points are in ascending order and are consecutive. So we can
simply subtract two encodings to get the number of representable values
between them (including one endpoint but not the other).
Unfortunately, the negative numbers are not adjacent and run the other
direction. To deal with this, if the number is negative, we transform
its encoding by subtracting from the encoding of -0. This gives us a
consecutive sequence of encodings from the greatest magnitude finite
negative number to the greatest finite number, in ascending order
except for wrapping at the maximum uint32_t value.
Note that this also maps the encoding of -0 to 0 (the encoding of +0),
so the two zeroes become one point, so they are counted only once.
*/
if (HighBit & ae) ae = HighBit - ae;
if (HighBit & be) be = HighBit - be;
// Return the distance between the two transformed encodings.
return be - ae;
}
static void Try(float a, float b)
{
printf("[%.99g, %.99g] contains %" PRIu32 " representable values.\n",
a, b, Distance(a, b) + 1);
}
int main(void)
{
if (sizeof(float) != sizeof(uint32_t))
{
fprintf(stderr, "Error, uint32_t must be the same size as float.\n");
exit(EXIT_FAILURE);
}
/* Prepare some test values: smallest positive (subnormal) value, largest
subnormal value, smallest normal value.
*/
float S1 = FLT_TRUE_MIN;
float N1 = FLT_MIN;
float S2 = N1 - S1;
// Test 0 <= a <= b.
Try( 0, 0);
Try( 0, S1);
Try( 0, S2);
Try( 0, N1);
Try( 0, 1./3);
Try(S1, S1);
Try(S1, S2);
Try(S1, N1);
Try(S1, 1./3);
Try(S2, S2);
Try(S2, N1);
Try(S2, 1./3);
Try(N1, N1);
Try(N1, 1./3);
// Test a <= b <= 0.
Try(-0., -0.);
Try(-S1, -0.);
Try(-S2, -0.);
Try(-N1, -0.);
Try(-1./3, -0.);
Try(-S1, -S1);
Try(-S2, -S1);
Try(-N1, -S1);
Try(-1./3, -S1);
Try(-S2, -S2);
Try(-N1, -S2);
Try(-1./3, -S2);
Try(-N1, -N1);
Try(-1./3, -N1);
// Test a <= 0 <= b.
Try(-0., +0.);
Try(-0., S1);
Try(-0., S2);
Try(-0., N1);
Try(-0., 1./3);
Try(-S1, +0.);
Try(-S1, S1);
Try(-S1, S2);
Try(-S1, N1);
Try(-S1, 1./3);
Try(-S2, +0.);
Try(-S2, S1);
Try(-S2, S2);
Try(-S2, N1);
Try(-S2, 1./3);
Try(-N1, +0.);
Try(-N1, S1);
Try(-N1, S2);
Try(-N1, N1);
Try(-1./3, 1./3);
Try(-1./3, +0.);
Try(-1./3, S1);
Try(-1./3, S2);
Try(-1./3, N1);
Try(-1./3, 1./3);
return 0;
}
