mirror of
https://github.com/ClickHouse/ClickHouse.git
synced 2024-11-23 16:12:01 +00:00
Double precision of geoDistance if the arguments are Float64
This commit is contained in:
Alexey Milovidov
parent
6323ca0960
commit
65650809b4
@ -896,6 +896,7 @@ class IColumn;
|
||||
M(Int64, ignore_cold_parts_seconds, 0, "Only available in ClickHouse Cloud. Exclude new data parts from SELECT queries until they're either pre-warmed (see cache_populated_by_fetch) or this many seconds old. Only for Replicated-/SharedMergeTree.", 0) \
|
||||
M(Int64, prefer_warmed_unmerged_parts_seconds, 0, "Only available in ClickHouse Cloud. If a merged part is less than this many seconds old and is not pre-warmed (see cache_populated_by_fetch), but all its source parts are available and pre-warmed, SELECT queries will read from those parts instead. Only for ReplicatedMergeTree. Note that this only checks whether CacheWarmer processed the part; if the part was fetched into cache by something else, it'll still be considered cold until CacheWarmer gets to it; if it was warmed, then evicted from cache, it'll still be considered warm.", 0) \
|
||||
M(Bool, iceberg_engine_ignore_schema_evolution, false, "Ignore schema evolution in Iceberg table engine and read all data using latest schema saved on table creation. Note that it can lead to incorrect result", 0) \
|
||||
M(Bool, geo_distance_returns_float64_on_float64_arguments, true, "If all four arguments to `geoDistance`, `greatCircleDistance`, `greatCircleAngle` functions are Float64, return Float64 and use double precision for internal calculations. In previous ClickHouse versions, the functions always returned Float32.", 0) \
|
||||
|
||||
// End of COMMON_SETTINGS
|
||||
// Please add settings related to formats into the FORMAT_FACTORY_SETTINGS, move obsolete settings to OBSOLETE_SETTINGS and obsolete format settings to OBSOLETE_FORMAT_SETTINGS.
|
||||
|
||||
@ -113,6 +113,7 @@ static std::map<ClickHouseVersion, SettingsChangesHistory::SettingsChanges> sett
|
||||
{"output_format_parquet_compression_method", "lz4", "zstd", "Parquet/ORC/Arrow support many compression methods, including lz4 and zstd. ClickHouse supports each and every compression method. Some inferior tools, such as 'duckdb', lack support for the faster `lz4` compression method, that's why we set zstd by default."},
|
||||
{"output_format_orc_compression_method", "lz4", "zstd", "Parquet/ORC/Arrow support many compression methods, including lz4 and zstd. ClickHouse supports each and every compression method. Some inferior tools, such as 'duckdb', lack support for the faster `lz4` compression method, that's why we set zstd by default."},
|
||||
{"output_format_pretty_highlight_digit_groups", false, true, "If enabled and if output is a terminal, highlight every digit corresponding to the number of thousands, millions, etc. with underline."},
|
||||
{"geo_distance_returns_float64_on_float64_arguments", false, true, "Increase the default precision."},
|
||||
}},
|
||||
{"24.2", {{"allow_suspicious_variant_types", true, false, "Don't allow creating Variant type with suspicious variants by default"},
|
||||
{"validate_experimental_and_suspicious_types_inside_nested_types", false, true, "Validate usage of experimental and suspicious types inside nested types"},
|
||||
|
||||
@ -42,121 +42,6 @@ namespace ErrorCodes
|
||||
namespace
|
||||
{
|
||||
|
||||
constexpr double PI = std::numbers::pi_v<double>;
|
||||
constexpr float PI_F = std::numbers::pi_v<float>;
|
||||
|
||||
constexpr float RAD_IN_DEG = static_cast<float>(PI / 180.0);
|
||||
constexpr float RAD_IN_DEG_HALF = static_cast<float>(PI / 360.0);
|
||||
|
||||
constexpr size_t COS_LUT_SIZE = 1024; // maxerr 0.00063%
|
||||
constexpr float COS_LUT_SIZE_F = 1024.0f; // maxerr 0.00063%
|
||||
constexpr size_t ASIN_SQRT_LUT_SIZE = 512;
|
||||
constexpr size_t METRIC_LUT_SIZE = 1024;
|
||||
|
||||
/** Earth radius in meters using WGS84 authalic radius.
|
||||
* We use this value to be consistent with H3 library.
|
||||
*/
|
||||
constexpr float EARTH_RADIUS = 6371007.180918475f;
|
||||
constexpr float EARTH_DIAMETER = 2 * EARTH_RADIUS;
|
||||
|
||||
|
||||
float cos_lut[COS_LUT_SIZE + 1]; /// cos(x) table
|
||||
float asin_sqrt_lut[ASIN_SQRT_LUT_SIZE + 1]; /// asin(sqrt(x)) * earth_diameter table
|
||||
|
||||
float sphere_metric_lut[METRIC_LUT_SIZE + 1]; /// sphere metric, unitless: the distance in degrees for one degree across longitude depending on latitude
|
||||
float sphere_metric_meters_lut[METRIC_LUT_SIZE + 1]; /// sphere metric: the distance in meters for one degree across longitude depending on latitude
|
||||
float wgs84_metric_meters_lut[2 * (METRIC_LUT_SIZE + 1)]; /// ellipsoid metric: the distance in meters across one degree latitude/longitude depending on latitude
|
||||
|
||||
|
||||
inline double sqr(double v)
|
||||
{
|
||||
return v * v;
|
||||
}
|
||||
|
||||
inline float sqrf(float v)
|
||||
{
|
||||
return v * v;
|
||||
}
|
||||
|
||||
void geodistInit()
|
||||
{
|
||||
for (size_t i = 0; i <= COS_LUT_SIZE; ++i)
|
||||
cos_lut[i] = static_cast<float>(cos(2 * PI * i / COS_LUT_SIZE)); // [0, 2 * pi] -> [0, COS_LUT_SIZE]
|
||||
|
||||
for (size_t i = 0; i <= ASIN_SQRT_LUT_SIZE; ++i)
|
||||
asin_sqrt_lut[i] = static_cast<float>(asin(
|
||||
sqrt(static_cast<double>(i) / ASIN_SQRT_LUT_SIZE))); // [0, 1] -> [0, ASIN_SQRT_LUT_SIZE]
|
||||
|
||||
for (size_t i = 0; i <= METRIC_LUT_SIZE; ++i)
|
||||
{
|
||||
double latitude = i * (PI / METRIC_LUT_SIZE) - PI * 0.5; // [-pi / 2, pi / 2] -> [0, METRIC_LUT_SIZE]
|
||||
|
||||
/// Squared metric coefficients (for the distance in meters) on a tangent plane, for latitude and longitude (in degrees),
|
||||
/// depending on the latitude (in radians).
|
||||
|
||||
/// https://github.com/mapbox/cheap-ruler/blob/master/index.js#L67
|
||||
wgs84_metric_meters_lut[i * 2] = static_cast<float>(sqr(111132.09 - 566.05 * cos(2 * latitude) + 1.20 * cos(4 * latitude)));
|
||||
wgs84_metric_meters_lut[i * 2 + 1] = static_cast<float>(sqr(111415.13 * cos(latitude) - 94.55 * cos(3 * latitude) + 0.12 * cos(5 * latitude)));
|
||||
|
||||
sphere_metric_meters_lut[i] = static_cast<float>(sqr((EARTH_DIAMETER * PI / 360) * cos(latitude)));
|
||||
|
||||
sphere_metric_lut[i] = static_cast<float>(sqr(cos(latitude)));
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
if (f > 180)
|
||||
f = 360 - f;
|
||||
return f;
|
||||
}
|
||||
|
||||
inline float geodistFastCos(float x)
|
||||
{
|
||||
float y = fabsf(x) * (COS_LUT_SIZE_F / PI_F / 2.0f);
|
||||
size_t i = floatToIndex(y);
|
||||
y -= i;
|
||||
i &= (COS_LUT_SIZE - 1);
|
||||
return cos_lut[i] + (cos_lut[i + 1] - cos_lut[i]) * y;
|
||||
}
|
||||
|
||||
inline float geodistFastSin(float x)
|
||||
{
|
||||
float y = fabsf(x) * (COS_LUT_SIZE_F / PI_F / 2.0f);
|
||||
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;
|
||||
}
|
||||
|
||||
/// fast implementation of asin(sqrt(x))
|
||||
/// max error in floats 0.00369%, in doubles 0.00072%
|
||||
inline float geodistFastAsinSqrt(float x)
|
||||
{
|
||||
if (x < 0.122f)
|
||||
{
|
||||
// distance under 4546 km, Taylor error under 0.00072%
|
||||
float y = sqrtf(x);
|
||||
return y + x * y * 0.166666666666666f + x * x * y * 0.075f + x * x * x * y * 0.044642857142857f;
|
||||
}
|
||||
if (x < 0.948f)
|
||||
{
|
||||
// distance under 17083 km, 512-entry LUT error under 0.00072%
|
||||
x *= ASIN_SQRT_LUT_SIZE;
|
||||
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
|
||||
}
|
||||
|
||||
|
||||
enum class Method
|
||||
{
|
||||
SPHERE_DEGREES,
|
||||
@ -164,18 +49,123 @@ enum class Method
|
||||
WGS84_METERS,
|
||||
};
|
||||
|
||||
}
|
||||
constexpr size_t ASIN_SQRT_LUT_SIZE = 512;
|
||||
constexpr size_t COS_LUT_SIZE = 1024; // maxerr 0.00063%
|
||||
constexpr size_t METRIC_LUT_SIZE = 1024;
|
||||
|
||||
template <typename T>
|
||||
struct Impl
|
||||
{
|
||||
static constexpr T PI = std::numbers::pi_v<T>;
|
||||
static constexpr T RAD_IN_DEG = static_cast<T>(PI / T(180.0));
|
||||
static constexpr T RAD_IN_DEG_HALF = static_cast<T>(PI / T(360.0));
|
||||
|
||||
static constexpr T COS_LUT_SIZE_F = T(1024.0);
|
||||
|
||||
/** Earth radius in meters using WGS84 authalic radius.
|
||||
* We use this value to be consistent with H3 library.
|
||||
*/
|
||||
static constexpr T EARTH_RADIUS = T(6371007.180918475);
|
||||
static constexpr T EARTH_DIAMETER = 2 * EARTH_RADIUS;
|
||||
|
||||
T cos_lut[COS_LUT_SIZE + 1]; /// cos(x) table
|
||||
T asin_sqrt_lut[ASIN_SQRT_LUT_SIZE + 1]; /// asin(sqrt(x)) * earth_diameter table
|
||||
|
||||
T sphere_metric_lut[METRIC_LUT_SIZE + 1]; /// sphere metric, unitless: the distance in degrees for one degree across longitude depending on latitude
|
||||
T sphere_metric_meters_lut[METRIC_LUT_SIZE + 1]; /// sphere metric: the distance in meters for one degree across longitude depending on latitude
|
||||
T wgs84_metric_meters_lut[2 * (METRIC_LUT_SIZE + 1)]; /// ellipsoid metric: the distance in meters across one degree latitude/longitude depending on latitude
|
||||
|
||||
static T sqr(T v) { return v * v; }
|
||||
|
||||
Impl()
|
||||
{
|
||||
for (size_t i = 0; i <= COS_LUT_SIZE; ++i)
|
||||
cos_lut[i] = std::cos(2 * PI * i / COS_LUT_SIZE); // [0, 2 * pi] -> [0, COS_LUT_SIZE]
|
||||
|
||||
for (size_t i = 0; i <= ASIN_SQRT_LUT_SIZE; ++i)
|
||||
asin_sqrt_lut[i] = std::asin(std::sqrt(static_cast<T>(i) / ASIN_SQRT_LUT_SIZE)); // [0, 1] -> [0, ASIN_SQRT_LUT_SIZE]
|
||||
|
||||
for (size_t i = 0; i <= METRIC_LUT_SIZE; ++i)
|
||||
{
|
||||
T latitude = i * (PI / METRIC_LUT_SIZE) - PI * T(0.5); // [-pi / 2, pi / 2] -> [0, METRIC_LUT_SIZE]
|
||||
|
||||
/// Squared metric coefficients (for the distance in meters) on a tangent plane, for latitude and longitude (in degrees),
|
||||
/// depending on the latitude (in radians).
|
||||
|
||||
/// https://github.com/mapbox/cheap-ruler/blob/master/index.js#L67
|
||||
wgs84_metric_meters_lut[i * 2] = sqr(T(111132.09) - T(566.05) * std::cos(T(2.0) * latitude) + T(1.20) * std::cos(T(4.0) * latitude));
|
||||
wgs84_metric_meters_lut[i * 2 + 1] = sqr(T(111415.13) * std::cos(latitude) - T(94.55) * std::cos(T(3.0) * latitude) + T(0.12) * std::cos(T(5.0) * latitude));
|
||||
sphere_metric_meters_lut[i] = sqr((EARTH_DIAMETER * PI / 360) * std::cos(latitude));
|
||||
|
||||
sphere_metric_lut[i] = sqr(std::cos(latitude));
|
||||
}
|
||||
}
|
||||
|
||||
static inline NO_SANITIZE_UNDEFINED size_t toIndex(T x)
|
||||
{
|
||||
/// Implementation specific behaviour on overflow or infinite value.
|
||||
return static_cast<size_t>(x);
|
||||
}
|
||||
|
||||
static inline T degDiff(T f)
|
||||
{
|
||||
f = std::abs(f);
|
||||
if (f > 180)
|
||||
f = 360 - f;
|
||||
return f;
|
||||
}
|
||||
|
||||
inline T fastCos(T x)
|
||||
{
|
||||
T y = std::abs(x) * (COS_LUT_SIZE_F / PI / T(2.0));
|
||||
size_t i = toIndex(y);
|
||||
y -= i;
|
||||
i &= (COS_LUT_SIZE - 1);
|
||||
return cos_lut[i] + (cos_lut[i + 1] - cos_lut[i]) * y;
|
||||
}
|
||||
|
||||
inline T fastSin(T x)
|
||||
{
|
||||
T y = std::abs(x) * (COS_LUT_SIZE_F / PI / T(2.0));
|
||||
size_t i = toIndex(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;
|
||||
}
|
||||
|
||||
/// fast implementation of asin(sqrt(x))
|
||||
/// max error in floats 0.00369%, in doubles 0.00072%
|
||||
inline T fastAsinSqrt(T x)
|
||||
{
|
||||
if (x < T(0.122))
|
||||
{
|
||||
// distance under 4546 km, Taylor error under 0.00072%
|
||||
T y = std::sqrt(x);
|
||||
return y + x * y * T(0.166666666666666) + x * x * y * T(0.075) + x * x * x * y * T(0.044642857142857);
|
||||
}
|
||||
if (x < T(0.948))
|
||||
{
|
||||
// distance under 17083 km, 512-entry LUT error under 0.00072%
|
||||
x *= ASIN_SQRT_LUT_SIZE;
|
||||
size_t i = toIndex(x);
|
||||
return asin_sqrt_lut[i] + (asin_sqrt_lut[i + 1] - asin_sqrt_lut[i]) * (x - i);
|
||||
}
|
||||
return std::asin(std::sqrt(x)); /// distance is over 17083 km, just compute exact
|
||||
}
|
||||
};
|
||||
|
||||
template <typename T> Impl<T> impl;
|
||||
|
||||
DECLARE_MULTITARGET_CODE(
|
||||
|
||||
namespace
|
||||
{
|
||||
|
||||
template <Method method>
|
||||
float distance(float lon1deg, float lat1deg, float lon2deg, float lat2deg)
|
||||
template <Method method, typename T>
|
||||
T distance(T lon1deg, T lat1deg, T lon2deg, T lat2deg)
|
||||
{
|
||||
float lat_diff = geodistDegDiff(lat1deg - lat2deg);
|
||||
float lon_diff = geodistDegDiff(lon1deg - lon2deg);
|
||||
T lat_diff = impl<T>.degDiff(lat1deg - lat2deg);
|
||||
T lon_diff = impl<T>.degDiff(lon1deg - lon2deg);
|
||||
|
||||
if (lon_diff < 13)
|
||||
{
|
||||
@ -187,51 +177,51 @@ float distance(float lon1deg, float lat1deg, float lon2deg, float lat2deg)
|
||||
/// (Remember how a plane flies from Amsterdam to New York)
|
||||
/// 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, METRIC_LUT_SIZE] indexes
|
||||
size_t latitude_midpoint_index = floatToIndex(latitude_midpoint) & (METRIC_LUT_SIZE - 1);
|
||||
T latitude_midpoint = (lat1deg + lat2deg + 180) * METRIC_LUT_SIZE / 360; // [-90, 90] degrees -> [0, METRIC_LUT_SIZE] indexes
|
||||
size_t latitude_midpoint_index = impl<T>.toIndex(latitude_midpoint) & (METRIC_LUT_SIZE - 1);
|
||||
|
||||
/// This is linear interpolation between two table items at index "latitude_midpoint_index" and "latitude_midpoint_index + 1".
|
||||
|
||||
float k_lat{};
|
||||
float k_lon{};
|
||||
T k_lat{};
|
||||
T k_lon{};
|
||||
|
||||
if constexpr (method == Method::SPHERE_DEGREES)
|
||||
{
|
||||
k_lat = 1;
|
||||
|
||||
k_lon = sphere_metric_lut[latitude_midpoint_index]
|
||||
+ (sphere_metric_lut[latitude_midpoint_index + 1] - sphere_metric_lut[latitude_midpoint_index]) * (latitude_midpoint - latitude_midpoint_index);
|
||||
k_lon = impl<T>.sphere_metric_lut[latitude_midpoint_index]
|
||||
+ (impl<T>.sphere_metric_lut[latitude_midpoint_index + 1] - impl<T>.sphere_metric_lut[latitude_midpoint_index]) * (latitude_midpoint - latitude_midpoint_index);
|
||||
}
|
||||
else if constexpr (method == Method::SPHERE_METERS)
|
||||
{
|
||||
k_lat = sqrf(EARTH_DIAMETER * PI_F / 360.0f);
|
||||
k_lat = impl<T>.sqr(impl<T>.EARTH_DIAMETER * impl<T>.PI / T(360.0));
|
||||
|
||||
k_lon = sphere_metric_meters_lut[latitude_midpoint_index]
|
||||
+ (sphere_metric_meters_lut[latitude_midpoint_index + 1] - sphere_metric_meters_lut[latitude_midpoint_index]) * (latitude_midpoint - latitude_midpoint_index);
|
||||
k_lon = impl<T>.sphere_metric_meters_lut[latitude_midpoint_index]
|
||||
+ (impl<T>.sphere_metric_meters_lut[latitude_midpoint_index + 1] - impl<T>.sphere_metric_meters_lut[latitude_midpoint_index]) * (latitude_midpoint - latitude_midpoint_index);
|
||||
}
|
||||
else if constexpr (method == Method::WGS84_METERS)
|
||||
{
|
||||
k_lat = wgs84_metric_meters_lut[latitude_midpoint_index * 2]
|
||||
+ (wgs84_metric_meters_lut[(latitude_midpoint_index + 1) * 2] - wgs84_metric_meters_lut[latitude_midpoint_index * 2]) * (latitude_midpoint - latitude_midpoint_index);
|
||||
k_lat = impl<T>.wgs84_metric_meters_lut[latitude_midpoint_index * 2]
|
||||
+ (impl<T>.wgs84_metric_meters_lut[(latitude_midpoint_index + 1) * 2] - impl<T>.wgs84_metric_meters_lut[latitude_midpoint_index * 2]) * (latitude_midpoint - latitude_midpoint_index);
|
||||
|
||||
k_lon = wgs84_metric_meters_lut[latitude_midpoint_index * 2 + 1]
|
||||
+ (wgs84_metric_meters_lut[(latitude_midpoint_index + 1) * 2 + 1] - wgs84_metric_meters_lut[latitude_midpoint_index * 2 + 1]) * (latitude_midpoint - latitude_midpoint_index);
|
||||
k_lon = impl<T>.wgs84_metric_meters_lut[latitude_midpoint_index * 2 + 1]
|
||||
+ (impl<T>.wgs84_metric_meters_lut[(latitude_midpoint_index + 1) * 2 + 1] - impl<T>.wgs84_metric_meters_lut[latitude_midpoint_index * 2 + 1]) * (latitude_midpoint - latitude_midpoint_index);
|
||||
}
|
||||
|
||||
/// Metric on a tangent plane: it differs from Euclidean metric only by scale of coordinates.
|
||||
return sqrtf(k_lat * lat_diff * lat_diff + k_lon * lon_diff * lon_diff);
|
||||
return std::sqrt(k_lat * lat_diff * lat_diff + k_lon * lon_diff * lon_diff);
|
||||
}
|
||||
else
|
||||
{
|
||||
// points too far away; use haversine
|
||||
|
||||
float a = sqrf(geodistFastSin(lat_diff * RAD_IN_DEG_HALF))
|
||||
+ geodistFastCos(lat1deg * RAD_IN_DEG) * geodistFastCos(lat2deg * RAD_IN_DEG) * sqrf(geodistFastSin(lon_diff * RAD_IN_DEG_HALF));
|
||||
T a = impl<T>.sqr(impl<T>.fastSin(lat_diff * impl<T>.RAD_IN_DEG_HALF))
|
||||
+ impl<T>.fastCos(lat1deg * impl<T>.RAD_IN_DEG) * impl<T>.fastCos(lat2deg * impl<T>.RAD_IN_DEG) * impl<T>.sqr(impl<T>.fastSin(lon_diff * impl<T>.RAD_IN_DEG_HALF));
|
||||
|
||||
if constexpr (method == Method::SPHERE_DEGREES)
|
||||
return (360.0f / PI_F) * geodistFastAsinSqrt(a);
|
||||
return (T(360.0) / impl<T>.PI) * impl<T>.fastAsinSqrt(a);
|
||||
else
|
||||
return EARTH_DIAMETER * geodistFastAsinSqrt(a);
|
||||
return impl<T>.EARTH_DIAMETER * impl<T>.fastAsinSqrt(a);
|
||||
}
|
||||
}
|
||||
|
||||
@ -241,13 +231,24 @@ template <Method method>
|
||||
class FunctionGeoDistance : public IFunction
|
||||
{
|
||||
public:
|
||||
static constexpr auto name =
|
||||
(method == Method::SPHERE_DEGREES) ? "greatCircleAngle"
|
||||
: ((method == Method::SPHERE_METERS) ? "greatCircleDistance"
|
||||
: "geoDistance");
|
||||
FunctionGeoDistance(ContextPtr context)
|
||||
{
|
||||
always_float32 = !context->getSettingsRef().geo_distance_returns_float64_on_float64_arguments;
|
||||
}
|
||||
|
||||
private:
|
||||
String getName() const override { return name; }
|
||||
bool always_float32;
|
||||
|
||||
String getName() const override
|
||||
{
|
||||
if constexpr (method == Method::SPHERE_DEGREES)
|
||||
return "greatCircleAngle";
|
||||
if constexpr (method == Method::SPHERE_METERS)
|
||||
return "greatCircleDistance";
|
||||
else
|
||||
return "geoDistance";
|
||||
}
|
||||
|
||||
size_t getNumberOfArguments() const override { return 4; }
|
||||
|
||||
bool useDefaultImplementationForConstants() const override { return true; }
|
||||
@ -255,22 +256,31 @@ private:
|
||||
|
||||
DataTypePtr getReturnTypeImpl(const DataTypes & arguments) const override
|
||||
{
|
||||
for (const auto arg_idx : collections::range(0, arguments.size()))
|
||||
bool has_float64 = false;
|
||||
|
||||
for (size_t arg_idx = 0; arg_idx < 4; ++arg_idx)
|
||||
{
|
||||
const auto * arg = arguments[arg_idx].get();
|
||||
if (!isNumber(WhichDataType(arg)))
|
||||
WhichDataType which(arguments[arg_idx]);
|
||||
|
||||
if (!isNumber(which))
|
||||
throw Exception(ErrorCodes::ILLEGAL_TYPE_OF_ARGUMENT, "Illegal type {} of argument {} of function {}. "
|
||||
"Must be numeric", arg->getName(), std::to_string(arg_idx + 1), getName());
|
||||
"Must be numeric", arguments[arg_idx]->getName(), std::to_string(arg_idx + 1), getName());
|
||||
|
||||
if (which.isFloat64())
|
||||
has_float64 = true;
|
||||
}
|
||||
|
||||
return std::make_shared<DataTypeFloat32>();
|
||||
if (has_float64 && !always_float32)
|
||||
return std::make_shared<DataTypeFloat64>();
|
||||
else
|
||||
return std::make_shared<DataTypeFloat32>();
|
||||
}
|
||||
|
||||
ColumnPtr executeImpl(const ColumnsWithTypeAndName & arguments, const DataTypePtr & result_type, size_t input_rows_count) const override
|
||||
{
|
||||
auto dst = ColumnVector<Float32>::create();
|
||||
auto & dst_data = dst->getData();
|
||||
dst_data.resize(input_rows_count);
|
||||
bool returns_float64 = WhichDataType(result_type).isFloat64();
|
||||
|
||||
auto dst = result_type->createColumn();
|
||||
|
||||
auto arguments_copy = arguments;
|
||||
for (auto & argument : arguments_copy)
|
||||
@ -280,10 +290,24 @@ private:
|
||||
argument.type = result_type;
|
||||
}
|
||||
|
||||
const auto * col_lon1 = convertArgumentColumnToFloat32(arguments_copy, 0);
|
||||
const auto * col_lat1 = convertArgumentColumnToFloat32(arguments_copy, 1);
|
||||
const auto * col_lon2 = convertArgumentColumnToFloat32(arguments_copy, 2);
|
||||
const auto * col_lat2 = convertArgumentColumnToFloat32(arguments_copy, 3);
|
||||
if (returns_float64)
|
||||
run<Float64>(arguments_copy, dst, input_rows_count);
|
||||
else
|
||||
run<Float32>(arguments_copy, dst, input_rows_count);
|
||||
|
||||
return dst;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void run(const ColumnsWithTypeAndName & arguments, MutableColumnPtr & dst, size_t input_rows_count) const
|
||||
{
|
||||
const auto * col_lon1 = convertArgumentColumn<T>(arguments, 0);
|
||||
const auto * col_lat1 = convertArgumentColumn<T>(arguments, 1);
|
||||
const auto * col_lon2 = convertArgumentColumn<T>(arguments, 2);
|
||||
const auto * col_lat2 = convertArgumentColumn<T>(arguments, 3);
|
||||
|
||||
auto & dst_data = assert_cast<ColumnVector<T> &>(*dst).getData();
|
||||
dst_data.resize(input_rows_count);
|
||||
|
||||
for (size_t row_num = 0; row_num < input_rows_count; ++row_num)
|
||||
{
|
||||
@ -291,20 +315,20 @@ private:
|
||||
col_lon1->getData()[row_num], col_lat1->getData()[row_num],
|
||||
col_lon2->getData()[row_num], col_lat2->getData()[row_num]);
|
||||
}
|
||||
|
||||
return dst;
|
||||
}
|
||||
|
||||
const ColumnFloat32 * convertArgumentColumnToFloat32(const ColumnsWithTypeAndName & arguments, size_t argument_index) const
|
||||
template <typename T>
|
||||
const ColumnVector<T> * convertArgumentColumn(const ColumnsWithTypeAndName & arguments, size_t argument_index) const
|
||||
{
|
||||
const auto * column_typed = checkAndGetColumn<ColumnFloat32>(arguments[argument_index].column.get());
|
||||
const auto * column_typed = checkAndGetColumn<ColumnVector<T>>(arguments[argument_index].column.get());
|
||||
if (!column_typed)
|
||||
throw Exception(
|
||||
ErrorCodes::ILLEGAL_COLUMN,
|
||||
"Illegal type {} of argument {} of function {}. Must be Float32.",
|
||||
"Illegal type {} of argument {} of function {}. Must be {}.",
|
||||
arguments[argument_index].type->getName(),
|
||||
argument_index + 1,
|
||||
getName());
|
||||
getName(),
|
||||
TypeName<T>);
|
||||
|
||||
return column_typed;
|
||||
}
|
||||
@ -316,18 +340,19 @@ template <Method method>
|
||||
class FunctionGeoDistance : public TargetSpecific::Default::FunctionGeoDistance<method>
|
||||
{
|
||||
public:
|
||||
explicit FunctionGeoDistance(ContextPtr context) : selector(context)
|
||||
explicit FunctionGeoDistance(ContextPtr context)
|
||||
: TargetSpecific::Default::FunctionGeoDistance<method>(context), selector(context)
|
||||
{
|
||||
selector.registerImplementation<TargetArch::Default,
|
||||
TargetSpecific::Default::FunctionGeoDistance<method>>();
|
||||
TargetSpecific::Default::FunctionGeoDistance<method>>(context);
|
||||
|
||||
#if USE_MULTITARGET_CODE
|
||||
selector.registerImplementation<TargetArch::AVX,
|
||||
TargetSpecific::AVX::FunctionGeoDistance<method>>();
|
||||
TargetSpecific::AVX::FunctionGeoDistance<method>>(context);
|
||||
selector.registerImplementation<TargetArch::AVX2,
|
||||
TargetSpecific::AVX2::FunctionGeoDistance<method>>();
|
||||
TargetSpecific::AVX2::FunctionGeoDistance<method>>(context);
|
||||
selector.registerImplementation<TargetArch::AVX512F,
|
||||
TargetSpecific::AVX512F::FunctionGeoDistance<method>>();
|
||||
TargetSpecific::AVX512F::FunctionGeoDistance<method>>(context);
|
||||
#endif
|
||||
}
|
||||
|
||||
@ -345,12 +370,13 @@ private:
|
||||
ImplementationSelector<IFunction> selector;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
REGISTER_FUNCTION(GeoDistance)
|
||||
{
|
||||
geodistInit();
|
||||
factory.registerFunction<FunctionGeoDistance<Method::SPHERE_DEGREES>>();
|
||||
factory.registerFunction<FunctionGeoDistance<Method::SPHERE_METERS>>();
|
||||
factory.registerFunction<FunctionGeoDistance<Method::WGS84_METERS>>();
|
||||
factory.registerFunction("greatCircleAngle", [](ContextPtr context) { return std::make_shared<FunctionGeoDistance<Method::SPHERE_DEGREES>>(std::move(context)); });
|
||||
factory.registerFunction("greatCircleDistance", [](ContextPtr context) { return std::make_shared<FunctionGeoDistance<Method::SPHERE_METERS>>(std::move(context)); });
|
||||
factory.registerFunction("geoDistance", [](ContextPtr context) { return std::make_shared<FunctionGeoDistance<Method::WGS84_METERS>>(std::move(context)); });
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -79,7 +79,7 @@ inline char * writeVarInt(Int64 x, char * ostr)
|
||||
return writeVarUInt(static_cast<UInt64>((x << 1) ^ (x >> 63)), ostr);
|
||||
}
|
||||
|
||||
namespace impl
|
||||
namespace varint_impl
|
||||
{
|
||||
|
||||
template <bool check_eof>
|
||||
@ -106,8 +106,8 @@ inline void readVarUInt(UInt64 & x, ReadBuffer & istr)
|
||||
inline void readVarUInt(UInt64 & x, ReadBuffer & istr)
|
||||
{
|
||||
if (istr.buffer().end() - istr.position() >= 10)
|
||||
return impl::readVarUInt<false>(x, istr);
|
||||
return impl::readVarUInt<true>(x, istr);
|
||||
return varint_impl::readVarUInt<false>(x, istr);
|
||||
return varint_impl::readVarUInt<true>(x, istr);
|
||||
}
|
||||
|
||||
inline void readVarUInt(UInt64 & x, std::istream & istr)
|
||||
|
||||
Loading…
Reference in New Issue
Block a user