mantaflow: Avoid computation errors by using hypot for vector norm…

CMakeLists.txt

@@ -725,22 +725,19 @@ endif()
 
 # if supported, use c++14 for all pre-processed c++ sources
 include(CheckCXXCompilerFlag)
+CHECK_CXX_COMPILER_FLAG("-std=c++23" COMPILER_SUPPORTS_CXX23)
 CHECK_CXX_COMPILER_FLAG("-std=c++14" COMPILER_SUPPORTS_CXX14)
 CHECK_CXX_COMPILER_FLAG("-std=c++11" COMPILER_SUPPORTS_CXX11)
 CHECK_CXX_COMPILER_FLAG("-std=c++0x" COMPILER_SUPPORTS_CXX0X)
 set(CXX_STANDARD_FLAG " ")
-if(COMPILER_SUPPORTS_CXX14)
-	set(CXX_STANDARD_FLAG " -std=c++14 ")
-elseif(COMPILER_SUPPORTS_CXX11)
-	set(CXX_STANDARD_FLAG " -std=c++11 ")
-elseif(COMPILER_SUPPORTS_CXX0X)
-	set(CXX_STANDARD_FLAG " -std=c++0x ")
+
+if(COMPILER_SUPPORTS_CXX23)
+	set(CXX_STANDARD_FLAG " -std=c++23 ")
 endif()
 set_source_files_properties(${PREPPED_SOURCES} ${PP_REGCPP} PROPERTIES COMPILE_FLAGS ${CXX_STANDARD_FLAG})
 
-# sanity check for openvdb support
-if(OPENVDB AND NOT COMPILER_SUPPORTS_CXX14)
-	message(FATAL_ERROR "Compilation with OpenVDB requires a compiler with C++14 support")
+if(NOT COMPILER_SUPPORTS_CXX23)
+	message(FATAL_ERROR "Hypot requires a compiler with C++23 support")
 endif()
 
 # necessary for newer LLVM clang versions, to prevent abundance of "Namify" template warnings... for compatibility disable unknown warnings

source/util/vector4d.h

@@ -260,6 +260,35 @@ inline bool operator!= (const Vector4D<S>& s1, const Vector4D<S>& s2) {
 // External functions
 //************************************************************************
 
+/* Based on libstdc++ implementation from:
+   https://github.com/gcc-mirror/gcc/blob/master/libstdc%2B%2B-v3/include/c_global/cmath#L3769
+*/
+template<typename _Tp>
+inline _Tp
+__hypot4( _Tp __x, _Tp __y, _Tp __z, _Tp __t ) {
+	__x = std::abs( __x );
+	__y = std::abs( __y );
+	__z = std::abs( __z );
+	__t = std::abs( __t );
+	if (_Tp __a = __x < __y ? (__y < __z ? (__z < __t ? __t : __z ) : (__y < __t ? __t : __y )) : __x < __z ? (__z < __t ? __t : __z ) : (__x < __t ? __t : __x ))
+		return __a * std::sqrt(( __x/__a ) * (__x/__a )
+								+ ( __y/__a ) * ( __y/__a )
+								+ ( __z/__a ) * ( __z/__a )
+								+ ( __t/__a ) * ( __t/__a ));
+	else
+		return {};
+}
+
+inline float
+hypot4( float __x, float __y, float __z, float __t ) {
+	return __hypot4<float>( __x, __y, __z, __t );
+}
+
+inline double
+hypot4( double __x, double __y, double __z, double __t ) {
+	return __hypot4<double>( __x, __y, __z, __t );
+}
+
 //! Dot product
 template<class S>
 inline S dot ( const Vector4D<S> &t, const Vector4D<S> &v ) {
@@ -280,8 +309,8 @@ inline Vector4D<S> cross ( const Vector4D<S> &t, const Vector4D<S> &v ) {
 //! Compute the magnitude (length) of the vector
 template<class S>
 inline S norm ( const Vector4D<S>& v ) {
-	S l = v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
-	return ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON ) ? 1. : sqrt ( l );
+	S l = hypot4( v.x, v.y, v.z, v.t );
+	return ( fabs( l - 1. ) < VECTOR_EPSILON ) ? 1. : l;
 }
 
 //! Compute squared magnitude
@@ -293,12 +322,11 @@ inline S normSquare ( const Vector4D<S>& v ) {
 //! Returns a normalized vector
 template<class S>
 inline Vector4D<S> getNormalized ( const Vector4D<S>& v ) {
-	S l = v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
-	if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON )
+	S l = hypot4( v.x, v.y, v.z, v.t );
+	if ( fabs( l - 1. ) < VECTOR_EPSILON )
 		return v; /* normalized "enough"... */
-	else if ( l > VECTOR_EPSILON*VECTOR_EPSILON )
-	{
-		S fac = 1./sqrt ( l );
+	else if ( l > VECTOR_EPSILON ) {
+		S fac = 1. / l;
 		return Vector4D<S> ( v.x*fac, v.y*fac, v.z*fac , v.t*fac );
 	}
 	else
@@ -310,12 +338,12 @@ inline Vector4D<S> getNormalized ( const Vector4D<S>& v ) {
 template<class S>
 inline S normalize ( Vector4D<S> &v ) {
 	S norm;
-	S l = v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
-	if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON ) {
+	S l = hypot4( v.x, v.y, v.z, v.t );
+	if ( fabs( l - 1. ) < VECTOR_EPSILON ) {
 		norm = 1.;
-	} else if ( l > VECTOR_EPSILON*VECTOR_EPSILON ) {
-		norm = sqrt ( l );
-		v *= 1./norm;        
+	} else if ( l > VECTOR_EPSILON ) {
+		norm = l;
+		v *= 1. / norm;
 	} else {
 		v = Vector4D<S>::Zero;
 		norm = 0.;

source/util/vectorbase.h

@@ -378,9 +378,9 @@ inline const Vector3D<S>& projectNormalTo ( const Vector3D<S>& v, const Vector3D
 //! (clamps to 0 and 1 with VECTOR_EPSILON)
 template<class S>
 inline S norm ( const Vector3D<S>& v ) {
-	S l = v.x*v.x + v.y*v.y + v.z*v.z;
-	if     (        l      <= VECTOR_EPSILON*VECTOR_EPSILON ) return(0.);
-	return ( fabs ( l-1. ) <  VECTOR_EPSILON*VECTOR_EPSILON ) ? 1. : sqrt ( l );
+	S l = std::hypot(v.x, v.y, v.z);
+	if     (        l      <= VECTOR_EPSILON) return (0.);
+	return ( fabs ( l-1. ) <  VECTOR_EPSILON) ? 1. : l;
 }
 
 //! Compute squared magnitude
@@ -398,12 +398,12 @@ inline Real normSquare(const int v) { return square(v); }
 //! Returns a normalized vector
 template<class S>
 inline Vector3D<S> getNormalized ( const Vector3D<S>& v ) {
-	S l = v.x*v.x + v.y*v.y + v.z*v.z;
-	if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON )
+	S l = std::hypot(v.x, v.y, v.z);
+	if ( fabs ( l-1. ) < VECTOR_EPSILON )
 		return v; /* normalized "enough"... */
-	else if ( l > VECTOR_EPSILON*VECTOR_EPSILON )
+	else if ( l > VECTOR_EPSILON )
 	{
-		S fac = 1./sqrt ( l );
+		S fac = 1./l;
 		return Vector3D<S> ( v.x*fac, v.y*fac, v.z*fac );
 	}
 	else
@@ -415,11 +415,11 @@ inline Vector3D<S> getNormalized ( const Vector3D<S>& v ) {
 template<class S>
 inline S normalize ( Vector3D<S> &v ) {
 	S norm;
-	S l = v.x*v.x + v.y*v.y + v.z*v.z;
-	if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON ) {
+	S l = std::hypot(v.x, v.y, v.z);
+	if ( fabs ( l-1. ) < VECTOR_EPSILON ) {
 		norm = 1.;
-	} else if ( l > VECTOR_EPSILON*VECTOR_EPSILON ) {
-		norm = sqrt ( l );
+	} else if ( l > VECTOR_EPSILON ) {
+		norm = l;
 		v *= 1./norm;        
 	} else {
 		v = Vector3D<S>::Zero;