Update vendored DuckDB sources to f8f3f4dcd95

CMakeLists.txt

@@ -152,6 +152,7 @@ set(DUCKDB_SRC_FILES
   src/duckdb/extension/parquet/parquet_file_metadata_cache.cpp
   src/duckdb/extension/parquet/parquet_float16.cpp
   src/duckdb/extension/parquet/parquet_geometry.cpp
+  src/duckdb/extension/parquet/parquet_int96.cpp
   src/duckdb/extension/parquet/parquet_metadata.cpp
   src/duckdb/extension/parquet/parquet_multi_file_info.cpp
   src/duckdb/extension/parquet/parquet_prefetch_cost_model.cpp

src/duckdb/extension/core_functions/aggregate/algebraic/avg.cpp

@@ -239,8 +239,9 @@ struct IntervalAverageOperation : public BaseSumOperation<AverageSetOperation, I
 			target.days += months_remainder / count;
 			days_remainder += months_remainder % count;
 
-			days_remainder *= Interval::MICROS_PER_DAY;
-			micros_remainder += days_remainder / count;
+			// the remainder in micros can exceed int64 for large counts, but divided by the count it is below two days
+			const auto shifted_days = hugeint_t(days_remainder) * hugeint_t(Interval::MICROS_PER_DAY);
+			micros_remainder += Hugeint::Cast<int64_t>(shifted_days / hugeint_t(count));
 			target.micros += micros_remainder;
 		}
 	}

src/duckdb/extension/core_functions/aggregate/distributive/arg_min_max.cpp

@@ -210,9 +210,7 @@ struct ArgMinMaxBase {
 		auto &context = input.GetClientContext();
 		auto &function = input.GetBoundFunction();
 		auto &arguments = input.GetArguments();
-		if (arguments[1]->GetReturnType().InternalType() == PhysicalType::VARCHAR) {
-			ExpressionBinder::PushCollation(context, arguments[1], arguments[1]->GetReturnType());
-		}
+		ExpressionBinder::PushCollation(context, arguments[1], arguments[1]->GetReturnType());
 		function.GetArguments()[0] = arguments[0]->GetReturnType();
 		function.GetArguments()[1] = arguments[1]->GetReturnType();
 		function.SetReturnType(arguments[0]->GetReturnType());
@@ -379,9 +377,7 @@ struct VectorArgMinMaxBase : ArgMinMaxBase<COMPARATOR> {
 		auto &context = input.GetClientContext();
 		auto &function = input.GetBoundFunction();
 		auto &arguments = input.GetArguments();
-		if (arguments[1]->GetReturnType().InternalType() == PhysicalType::VARCHAR) {
-			ExpressionBinder::PushCollation(context, arguments[1], arguments[1]->GetReturnType());
-		}
+		ExpressionBinder::PushCollation(context, arguments[1], arguments[1]->GetReturnType());
 		function.GetArguments()[0] = arguments[0]->GetReturnType();
 		function.GetArguments()[1] = arguments[1]->GetReturnType();
 		function.SetReturnType(arguments[0]->GetReturnType());
@@ -886,13 +882,15 @@ void SpecializeArgMinMaxNullNFunction(PhysicalType val_type, PhysicalType arg_ty
 
 template <ArgMinMaxNullHandling NULL_HANDLING, bool NULLS_LAST, class COMPARATOR>
 unique_ptr<FunctionData> ArgMinMaxNBind(BindAggregateFunctionInput &input) {
+	auto &context = input.GetClientContext();
 	auto &function = input.GetBoundFunction();
 	auto &arguments = input.GetArguments();
 	for (auto &arg : arguments) {
 		if (arg->GetReturnType().id() == LogicalTypeId::UNKNOWN) {
 			throw ParameterNotResolvedException();
 		}
 	}
+	ExpressionBinder::PushCollation(context, arguments[1], arguments[1]->GetReturnType());
 
 	const auto val_type = arguments[0]->GetReturnType().InternalType();
 	const auto arg_type = arguments[1]->GetReturnType().InternalType();

src/duckdb/extension/core_functions/aggregate/distributive/bitstring_agg.cpp

@@ -123,15 +123,12 @@ struct BitStringAggOperation {
 		if (min > max) {
 			throw InvalidInputException("Invalid explicit bitstring range: Minimum (%d) > maximum (%d)", min, max);
 		}
-		INPUT_TYPE result;
-		if (!TrySubtractOperator::Operation(max, min, result)) {
+		// the difference of two (at most 64-bit) values does not always fit in their own type, but it fits in a hugeint
+		const auto range = Hugeint::Convert(max) - Hugeint::Convert(min);
+		if (range >= Hugeint::Convert(NumericLimits<idx_t>::Maximum())) {
 			return NumericLimits<idx_t>::Maximum();
 		}
-		auto val = NumericCast<idx_t>(result);
-		if (val == NumericLimits<idx_t>::Maximum()) {
-			return val;
-		}
-		return val + 1;
+		return Hugeint::Cast<idx_t>(range) + 1;
 	}
 
 	template <class INPUT_TYPE, class STATE>

src/duckdb/extension/core_functions/aggregate/distributive/kurtosis.cpp

@@ -1,3 +1,4 @@
+#include "duckdb/common/operator/add.hpp"
 #include "core_functions/aggregate/distributive_functions.hpp"
 #include "duckdb/common/exception.hpp"
 #include "duckdb/common/vector_operations/vector_operations.hpp"
@@ -8,15 +9,17 @@ namespace duckdb {
 
 namespace {
 
+//! Online central moments (Terriberry / Pébay, DOI: 10.2172/1028931), the higher-moment extension of Welford's
+//! method used by var_samp: M2..M4 are exactly zero for constant input instead of rounding noise of power sums
 struct KurtosisState {
-	static constexpr const char *STATE_NAMES[] = {"n", "sum", "sum_sqr", "sum_cub", "sum_four"};
+	static constexpr const char *STATE_NAMES[] = {"n", "mean", "m2", "m3", "m4"};
 	using STATE_TYPE = StructStateType<idx_t, double, double, double, double>;
 
 	idx_t n;
-	double sum;
-	double sum_sqr;
-	double sum_cub;
-	double sum_four;
+	double mean; //  M1
+	double m2;   //  M2 = sum((x - mean)^2)
+	double m3;   //  M3 = sum((x - mean)^3)
+	double m4;   //  M4 = sum((x - mean)^4)
 };
 
 struct KurtosisFlagBiasCorrection {};
@@ -28,64 +31,82 @@ struct KurtosisOperation {
 	template <class INPUT_TYPE, class STATE, class OP>
 	static void ConstantOperation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &unary_input,
 	                              idx_t count) {
-		for (idx_t i = 0; i < count; i++) {
-			Operation<INPUT_TYPE, STATE, OP>(state, input, unary_input);
-		}
+		// a block of identical values has mean "input" and zero central moments
+		STATE constant;
+		constant.n = count;
+		constant.mean = static_cast<double>(input);
+		constant.m2 = 0;
+		constant.m3 = 0;
+		constant.m4 = 0;
+		Combine<STATE, OP>(constant, state, unary_input.input);
 	}
 
 	template <class INPUT_TYPE, class STATE, class OP>
-	static void Operation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &unary_input) {
+	static void Operation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &) {
+		const double n1 = static_cast<double>(state.n);
 		state.n++;
-		state.sum += input;
-		state.sum_sqr += pow(input, 2);
-		state.sum_cub += pow(input, 3);
-		state.sum_four += pow(input, 4);
+		const double n = static_cast<double>(state.n);
+		const double delta = static_cast<double>(input) - state.mean;
+		const double delta_n = delta / n;
+		const double delta_n2 = delta_n * delta_n;
+		const double term1 = delta * delta_n * n1;
+		state.mean += delta_n;
+		state.m4 += term1 * delta_n2 * (n * n - 3 * n + 3) + 6 * delta_n2 * state.m2 - 4 * delta_n * state.m3;
+		state.m3 += term1 * delta_n * (n - 2) - 3 * delta_n * state.m2;
+		state.m2 += term1;
 	}
 
 	template <class STATE, class OP>
 	static void Combine(const STATE &source, STATE &target, AggregateInputData &) {
 		if (source.n == 0) {
 			return;
 		}
-		target.n += source.n;
-		target.sum += source.sum;
-		target.sum_sqr += source.sum_sqr;
-		target.sum_cub += source.sum_cub;
-		target.sum_four += source.sum_four;
+		if (target.n == 0) {
+			target = source;
+			return;
+		}
+		idx_t count;
+		if (!TryAddOperator::Operation(target.n, source.n, count)) {
+			throw OutOfRangeException("Overflow in aggregate state combine");
+		}
+		const double na = static_cast<double>(target.n);
+		const double nb = static_cast<double>(source.n);
+		const double n = static_cast<double>(count);
+		const double delta = source.mean - target.mean;
+		const double delta_n = delta / n;
+		const double delta_n2 = delta_n * delta_n;
+		target.m4 += source.m4 + delta * delta_n * delta_n2 * na * nb * (na * na - na * nb + nb * nb) +
+		             6 * delta_n2 * (na * na * source.m2 + nb * nb * target.m2) +
+		             4 * delta_n * (na * source.m3 - nb * target.m3);
+		target.m3 +=
+		    source.m3 + delta * delta_n2 * na * nb * (na - nb) + 3 * delta_n * (na * source.m2 - nb * target.m2);
+		target.m2 += source.m2 + delta * delta_n * na * nb;
+		target.mean += delta_n * nb;
+		target.n = count;
 	}
 
 	template <class TARGET_TYPE, class STATE>
 	static void Finalize(STATE &state, TARGET_TYPE &target, AggregateFinalizeData &finalize_data) {
-		auto n = (double)state.n;
-		if (n <= 1) {
+		if (state.n <= 1) {
 			finalize_data.ReturnNull();
 			return;
 		}
-		if (std::is_same<KURTOSIS_FLAG, KurtosisFlagBiasCorrection>::value && n <= 3) {
+		if (std::is_same<KURTOSIS_FLAG, KurtosisFlagBiasCorrection>::value && state.n <= 3) {
 			finalize_data.ReturnNull();
 			return;
 		}
-		double temp = 1 / n;
-		//! This is necessary due to linux 32 bits
-		long double temp_aux = 1 / n;
-		if (state.sum_sqr - state.sum * state.sum * temp == 0 ||
-		    state.sum_sqr - state.sum * state.sum * temp_aux == 0) {
-			finalize_data.ReturnNull();
-			return;
-		}
-		double m4 =
-		    temp * (state.sum_four - 4 * state.sum_cub * state.sum * temp +
-		            6 * state.sum_sqr * state.sum * state.sum * temp * temp - 3 * pow(state.sum, 4) * pow(temp, 3));
-
-		double m2 = temp * (state.sum_sqr - state.sum * state.sum * temp);
-		if (m2 <= 0) { // m2 shouldn't be below 0 but floating points are weird
+		// kurtosis is undefined for zero variance (M2 is exactly 0 for constant input)
+		if (state.m2 <= 0) {
 			finalize_data.ReturnNull();
 			return;
 		}
+		const double n = static_cast<double>(state.n);
+		// m4 / m2^2 with the population moments m4 = M4 / n and m2 = M2 / n
+		const double m4_over_m2_sqr = n * (state.m4 / state.m2) / state.m2;
 		if (std::is_same<KURTOSIS_FLAG, KurtosisFlagNoBiasCorrection>::value) {
-			target = m4 / (m2 * m2) - 3;
+			target = m4_over_m2_sqr - 3;
 		} else {
-			target = (n - 1) * ((n + 1) * m4 / (m2 * m2) - 3 * (n - 1)) / ((n - 2) * (n - 3));
+			target = (n - 1) * ((n + 1) * m4_over_m2_sqr - 3 * (n - 1)) / ((n - 2) * (n - 3));
 		}
 		if (!Value::DoubleIsFinite(target)) {
 			throw OutOfRangeException("Kurtosis is out of range!");

src/duckdb/extension/core_functions/aggregate/distributive/skew.cpp

@@ -3,49 +3,70 @@
 #include "duckdb/common/vector_operations/vector_operations.hpp"
 #include "duckdb/planner/expression/bound_aggregate_expression.hpp"
 #include "duckdb/common/algorithm.hpp"
-#include <limits>
 
 namespace duckdb {
 
 namespace {
 
+//! Online central moments (Terriberry / Pébay, DOI: 10.2172/1028931), the higher-moment extension of Welford's
+//! method used by var_samp: M2 and M3 are exactly zero for constant input instead of rounding noise of power sums
 struct SkewState {
-	static constexpr const char *STATE_NAMES[] = {"n", "sum", "sum_sqr", "sum_cub"};
+	static constexpr const char *STATE_NAMES[] = {"n", "mean", "m2", "m3"};
 	using STATE_TYPE = StructStateType<idx_t, double, double, double>;
 
 	idx_t n;
-	double sum;
-	double sum_sqr;
-	double sum_cub;
+	double mean; //  M1
+	double m2;   //  M2 = sum((x - mean)^2)
+	double m3;   //  M3 = sum((x - mean)^3)
 };
 
 struct SkewnessOperation {
 	template <class INPUT_TYPE, class STATE, class OP>
 	static void ConstantOperation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &unary_input,
 	                              idx_t count) {
-		for (idx_t i = 0; i < count; i++) {
-			Operation<INPUT_TYPE, STATE, OP>(state, input, unary_input);
-		}
+		// a block of identical values has mean "input" and zero central moments
+		STATE constant;
+		constant.n = count;
+		constant.mean = static_cast<double>(input);
+		constant.m2 = 0;
+		constant.m3 = 0;
+		Combine<STATE, OP>(constant, state, unary_input.input);
 	}
 
 	template <class INPUT_TYPE, class STATE, class OP>
-	static void Operation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &unary_input) {
+	static void Operation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &) {
+		const double n1 = static_cast<double>(state.n);
 		state.n++;
-		state.sum += input;
-		state.sum_sqr += pow(input, 2);
-		state.sum_cub += pow(input, 3);
+		const double n = static_cast<double>(state.n);
+		const double delta = static_cast<double>(input) - state.mean;
+		const double delta_n = delta / n;
+		const double term1 = delta * delta_n * n1;
+		state.mean += delta_n;
+		state.m3 += term1 * delta_n * (n - 2) - 3 * delta_n * state.m2;
+		state.m2 += term1;
 	}
 
 	template <class STATE, class OP>
 	static void Combine(const STATE &source, STATE &target, AggregateInputData &) {
 		if (source.n == 0) {
 			return;
 		}
-
-		target.n += source.n;
-		target.sum += source.sum;
-		target.sum_sqr += source.sum_sqr;
-		target.sum_cub += source.sum_cub;
+		if (target.n == 0) {
+			target = source;
+			return;
+		}
+		const auto count = target.n + source.n;
+		D_ASSERT(count >= target.n); // This is a check that we are not overflowing
+		const double na = static_cast<double>(target.n);
+		const double nb = static_cast<double>(source.n);
+		const double n = static_cast<double>(count);
+		const double delta = source.mean - target.mean;
+		const double delta_n = delta / n;
+		target.m3 += source.m3 + delta * delta_n * delta_n * na * nb * (na - nb) +
+		             3 * delta_n * (na * source.m2 - nb * target.m2);
+		target.m2 += source.m2 + delta * delta_n * na * nb;
+		target.mean += delta_n * nb;
+		target.n = count;
 	}
 
 	template <class TARGET_TYPE, class STATE>
@@ -54,26 +75,15 @@ struct SkewnessOperation {
 			finalize_data.ReturnNull();
 			return;
 		}
-		double n = state.n;
-		double temp = 1 / n;
-		double raw_m2 = state.sum_sqr - state.sum * state.sum * temp;
-		// Only treat finite second-moment noise as zero; overflow should still surface as out-of-range.
-		if (Value::DoubleIsFinite(raw_m2) && Value::DoubleIsFinite(state.sum_sqr) &&
-		    // Scale the tolerance by the accumulated squared magnitude instead of using a fixed epsilon.
-		    std::abs(raw_m2) <= std::numeric_limits<double>::epsilon() * std::max(1.0, std::abs(state.sum_sqr))) {
-			finalize_data.ReturnNull();
-			return;
-		}
-		double variance = temp * raw_m2;
-		if (variance <= 0) {
+		// skewness is undefined for zero variance (M2 is exactly 0 for constant input)
+		if (state.m2 <= 0) {
 			finalize_data.ReturnNull();
 			return;
 		}
-		auto p = std::pow(variance, 3);
-		double div = std::sqrt(p);
-		double temp1 = std::sqrt(n * (n - 1)) / (n - 2);
-		target = temp1 * temp *
-		         (state.sum_cub - 3 * state.sum_sqr * state.sum * temp + 2 * pow(state.sum, 3) * temp * temp) / div;
+		const double n = static_cast<double>(state.n);
+		// g1 = (M3 / n) / (M2 / n)^1.5 = sqrt(n) * M3 / M2^1.5, with the bias correction sqrt(n (n - 1)) / (n - 2)
+		const double g1 = std::sqrt(n) * (state.m3 / state.m2) / std::sqrt(state.m2);
+		target = std::sqrt(n * (n - 1)) / (n - 2) * g1;
 		if (!Value::DoubleIsFinite(target)) {
 			throw OutOfRangeException("SKEW is out of range!");
 		}