Additions

This commit is contained in:
Alexey Milovidov 2021-01-23 03:14:16 +03:00
parent 0528ca60d6
commit 9a28823041

View File

@ -99,6 +99,12 @@ void geodistInit()
}
}
inline NO_SANITIZE_UNDEFINED size_t floatToIndex(float x)
{
/// Implementation specific behaviour on overflow or infinite value.
return static_cast<size_t>(x);
}
inline float geodistDegDiff(float f)
{
f = fabsf(f);
@ -110,7 +116,7 @@ inline float geodistDegDiff(float f)
inline float geodistFastCos(float x)
{
float y = fabsf(x) * (COS_LUT_SIZE / PI / 2);
size_t i = static_cast<size_t>(y);
size_t i = floatToIndex(y);
y -= i;
i &= (COS_LUT_SIZE - 1);
return cos_lut[i] + (cos_lut[i + 1] - cos_lut[i]) * y;
@ -119,7 +125,7 @@ inline float geodistFastCos(float x)
inline float geodistFastSin(float x)
{
float y = fabsf(x) * (COS_LUT_SIZE / PI / 2);
size_t i = static_cast<size_t>(y);
size_t i = floatToIndex(y);
y -= i;
i = (i - COS_LUT_SIZE / 4) & (COS_LUT_SIZE - 1); // cos(x - pi / 2) = sin(x), costable / 4 = pi / 2
return cos_lut[i] + (cos_lut[i + 1] - cos_lut[i]) * y;
@ -139,7 +145,7 @@ inline float geodistFastAsinSqrt(float x)
{
// distance under 17083 km, 512-entry LUT error under 0.00072%
x *= ASIN_SQRT_LUT_SIZE;
size_t i = static_cast<size_t>(x);
size_t i = floatToIndex(x);
return asin_sqrt_lut[i] + (asin_sqrt_lut[i + 1] - asin_sqrt_lut[i]) * (x - i);
}
return asinf(sqrtf(x)); // distance over 17083 km, just compute exact
@ -160,12 +166,6 @@ DECLARE_MULTITARGET_CODE(
namespace
{
inline NO_SANITIZE_UNDEFINED size_t metricLUTIndex(float latitude_midpoint)
{
/// Implementation specific behaviour on overflow or infinite value.
return static_cast<size_t>(latitude_midpoint) & (METRIC_LUT_SIZE - 1);
}
template <Method method>
float distance(float lon1deg, float lat1deg, float lon2deg, float lat2deg)
{
@ -183,7 +183,7 @@ float distance(float lon1deg, float lat1deg, float lon2deg, float lat2deg)
/// But if longitude is close but latitude is different enough, there is no difference between meridian and great circle line.
float latitude_midpoint = (lat1deg + lat2deg + 180) * METRIC_LUT_SIZE / 360; // [-90, 90] degrees -> [0, KTABLE] indexes
size_t latitude_midpoint_index = metricLUTIndex(latitude_midpoint);
size_t latitude_midpoint_index = floatToIndex(latitude_midpoint) & (METRIC_LUT_SIZE - 1);
/// This is linear interpolation between two table items at index "latitude_midpoint_index" and "latitude_midpoint_index + 1".