geohash_encode() and geohash_fast_encode() produces different hashes in some cases

[Happlo]

Here are two test cases where the function gehoash_encode() and geohash_fast_encode() produces different result. In the last case geohash_fast_encode() is using too many bits.

GeoHashRange const LAT_RANGE{ 90.0,-90.0 };
GeoHashRange const LON_RANGE{ 180.0, -180.0 };
GeoHashBits hash_0, hash_0_fast, hash_180, hash_180_fast;
double lat_point = 75.0; 
double lon_point_0 = 0.0;
double lon_point_180 = 180.0;
uint8_t steps = 4;

geohash_encode(LAT_RANGE, LON_RANGE, lat_point, lon_point_0, steps, &hash_0);
geohash_fast_encode(LAT_RANGE, LON_RANGE, lat_point, lon_point_0, steps, &hash_0_fast);
// I expect hash_0.bits (0x7e) and hash_0_fast.bits (0xd4) to be equal

geohash_encode(LAT_RANGE, LON_RANGE, lat_point, lon_point_180, steps, &hash_180);
geohash_fast_encode(LAT_RANGE, LON_RANGE, lat_point, lon_point_180, steps, &hash_180_fast);
// I expect hash_180.bits (0xfe) and hash_180_fast.bits (0x254) to be equal
// I expect hash_180_fast.bits to be using 8 bits and not 10.

I think the problem is when mapping position from [range.min, range.max] to [0, 2^steps]. If the position is mapped to a whole number it should be decremented by one to get same behaviour as geohash_encode(). For example, if step is 4 and offset is 16.0 (maximum offset) we use 15 instead of 16 since 16 requires 5 bits. I think it is fixed with the following code, however I am not sure about the performance impact. A more efficient solution is probably possible.

static int is_whole(double number)
{
    return floor(number) == number;
}

static uint32_t calculate_offset(double position, GeoHashRange range, uint8_t step)
{
    // map position from [range.min, range.max] to [0, 2^steps]
    double offset = (position - range.min) * (1LL << step) / (range.max - range.min);
    if (offset > 0.0 && is_whole(offset))
    {
        // If the offset is a whole number, decrement it by one to get same behavior as the normal
        // gehoash_encode(). This is especially important if for example step is 4 and offset is 16.0
        // (maximum offset). We cannot use 16 when step is 4 since it is 5 bits.
        return ((uint32_t) offset) - 1;
    }
    return (uint32_t)offset;
}

int geohash_fast_encode(
        GeoHashRange lat_range, GeoHashRange lon_range, double latitude,
        double longitude, uint8_t step,  GeoHashBits* hash)
{
    if (NULL == hash || step > 32 || step == 0)
    {
        return -1;
    }
    hash->bits = 0;
    hash->step = step;
    if   (latitude < lat_range.min || latitude > lat_range.max
      || longitude < lon_range.min || longitude > lon_range.max)
    {
        return -1;
    }

    // The algorithm computes the morton code for the geohash location within
    // the range this can be done MUCH more efficiently using the following code

    uint32_t ilato = calculate_offset(latitude, lat_range, step);
    uint32_t ilono = calculate_offset(longitude, lon_range, step);

    //interleave the bits to create the morton code.  No branching and no bounding
    hash->bits = interleave64(ilato, ilono);
    return 0;
}