diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Do_intersect_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Do_intersect_3.h index 960129a237b3..2062c3ddecee 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Do_intersect_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Do_intersect_3.h @@ -345,7 +345,7 @@ class Do_intersect_3 double aprz = CGAL::abs(prz); double apsz = CGAL::abs(psz); -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX CGAL::Max mmax; maxx = mmax(maxx, aprx, apsx); @@ -365,14 +365,14 @@ class Do_intersect_3 double eps = 5.1107127829973299e-15 * maxx * maxy * maxz; -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX #if 0 CGAL::Min mmin; double tmp = mmin(maxx, maxy, maxz); maxz = mmax(maxx, maxy, maxz); maxx = tmp; #else - sse2minmax(maxx,maxy,maxz); + simd_minmax(maxx,maxy,maxz); // maxy can contain ANY element #endif #else diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Orientation_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Orientation_3.h index ce73ed3e83c5..38f90e3353fe 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Orientation_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Orientation_3.h @@ -85,7 +85,7 @@ class Orientation_3 double aprz = CGAL::abs(prz); double apsz = CGAL::abs(psz); -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX CGAL::Max mmax; maxx = mmax(maxx, aprx, apsx); @@ -105,14 +105,14 @@ class Orientation_3 double eps = 5.1107127829973299e-15 * maxx * maxy * maxz; -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX #if 0 CGAL::Min mmin; double tmp = mmin(maxx, maxy, maxz); maxz = mmax(maxx, maxy, maxz); maxx = tmp; #else - sse2minmax(maxx,maxy,maxz); + simd_minmax(maxx,maxy,maxz); // maxy can contain ANY element #endif #else diff --git a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h index 90b128c5eee6..a93940fc4743 100644 --- a/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h +++ b/Filtered_kernel/include/CGAL/Filtered_kernel/internal/Static_filters/Side_of_oriented_sphere_3.h @@ -86,7 +86,7 @@ class Side_of_oriented_sphere_3 double artz = CGAL::abs(rtz); double astz = CGAL::abs(stz); -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX CGAL::Max mmax; maxx = mmax(maxx, aqtx, artx, astx); maxy = mmax(maxy, aqty, arty, asty); @@ -107,14 +107,14 @@ class Side_of_oriented_sphere_3 double eps = 1.2466136531027298e-13 * maxx * maxy * maxz; -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX /* CGAL::Min mmin; double tmp = mmin(maxx, maxy, maxz); maxz = mmax(maxx, maxy, maxz); maxx = tmp; */ - sse2minmax(maxx,maxy,maxz); + simd_minmax(maxx,maxy,maxz); // maxy can contain ANY element #else diff --git a/Mesh_3/include/CGAL/Mesh_3/Robust_intersection_traits_3.h b/Mesh_3/include/CGAL/Mesh_3/Robust_intersection_traits_3.h index f4eeebecef89..13f10073eb39 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Robust_intersection_traits_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Robust_intersection_traits_3.h @@ -120,7 +120,7 @@ struct Vector_plane_orientation_3_static_filter : double aprz = CGAL::abs(acz); double apsz = CGAL::abs(pqz); -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX CGAL::Max mmax; maxx = mmax(maxx, aprx, apsx); @@ -140,14 +140,14 @@ struct Vector_plane_orientation_3_static_filter : double eps = 5.1107127829973299e-15 * maxx * maxy * maxz; -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX #if 0 CGAL::Min mmin; double tmp = mmin(maxx, maxy, maxz); maxz = mmax(maxx, maxy, maxz); maxx = tmp; #else - sse2minmax(maxx,maxy,maxz); + simd_minmax(maxx,maxy,maxz); // maxy can contain ANY element #endif #else diff --git a/Number_types/include/CGAL/FPU.h b/Number_types/include/CGAL/FPU.h index 13b02adec19f..422c2d0f25b4 100644 --- a/Number_types/include/CGAL/FPU.h +++ b/Number_types/include/CGAL/FPU.h @@ -44,6 +44,9 @@ extern "C" { # else # include # endif +#elif defined _MSC_VER && defined _M_ARM64 + // MSVC on ARM64: _controlfp_s is x86-only; fenv.h is the correct interface. +# include #elif defined _MSC_VER || defined __sparc__ || \ (defined __i386__ && !defined __PGI && !defined __SUNPRO_CC \ && !defined __SSE2__) @@ -119,6 +122,15 @@ extern "C" { !defined CGAL_ALWAYS_ROUND_TO_NEAREST # define CGAL_USE_SSE2 1 #endif + +#if defined(__aarch64__) || defined(_M_ARM64) +# include +# define CGAL_HAS_NEON 1 +#endif + +#if defined(CGAL_HAS_NEON) && !defined(CGAL_ALWAYS_ROUND_TO_NEAREST) +# define CGAL_USE_NEON 1 +#endif #ifdef CGAL_CFG_DENORMALS_COMPILE_BUG double& get_static_minimin(); // Defined in Interval_arithmetic_impl.h #endif @@ -283,7 +295,49 @@ inline __m128d swap_m128d(__m128d x){ return _mm_shuffle_pd(x, x, 1); # endif } +#endif // CGAL_USE_SSE2 + +// NEON opacify & swap helpers +#ifdef CGAL_USE_NEON + +// Prevent the compiler from reordering or constant-folding across a +// fesetround() boundary, same role as IA_opacify128 on SSE2. +inline float64x2_t IA_opacify_neon(float64x2_t x) +{ +#if defined __GNUG__ || defined __clang__ + // "+w" pins the value in a NEON/FP register. + asm volatile ("" : "+w"(x)); + return x; +#else + // MSVC on ARM64: no inline asm available. + // We use the same volatile+memcpy. + volatile float64x2_t e = x; + std::memcpy(&x, (void*)&e, 16); + return x; +#endif +} + +// Weaker barrier: blocks algebraic reassociation but not fesetround migration. +inline float64x2_t IA_opacify_neon_weak(float64x2_t x) +{ +#if defined __GNUG__ || defined __clang__ + asm ("" : "+mw"(x)); + return x; +#else + // MSVC: fall back to the strong barrier; no weaker option without asm. + return IA_opacify_neon(x); #endif +} + +// Swap the two double lanes: {a, b} -> {b, a} +// Counterpart of swap_m128d(). vextq_f64 with offset 1 rotates the two +// 64-bit elements, which on a 2-element vector is equivalent to a swap. +inline float64x2_t swap_f64x2(float64x2_t x) +{ + return vextq_f64(x, x, 1); +} + +#endif // CGAL_USE_NEON // Interval arithmetic needs to protect against double-rounding effects // caused by excess FPU precision, even if it forces the 53bit mantissa @@ -365,7 +419,20 @@ inline double CGAL_IA_SQUARE(double a){ double b = CGAL_IA_STOP_CPROP(a); // only once return IA_up(b*b); } -#if defined CGAL_SAFE_SSE2 + +// Handle MSVC ARM64 first because SSE2 intrinsics and _controlfp_s are x86-only. +// This prevents ARM64 from incorrectly selecting unsupported x86 code paths. +#if defined ( _MSC_VER ) && defined ( _M_ARM64 ) +// MSVC on ARM64: is already included at the top of this file. +#define CGAL_IA_SETFPCW(CW) fesetround(CW) +#define CGAL_IA_GETFPCW(CW) CW = fegetround() +typedef int FPU_CW_t; +#define CGAL_FE_TONEAREST FE_TONEAREST +#define CGAL_FE_TOWARDZERO FE_TOWARDZERO +#define CGAL_FE_UPWARD FE_UPWARD +#define CGAL_FE_DOWNWARD FE_DOWNWARD + +#elif defined CGAL_SAFE_SSE2 #define CGAL_IA_SETFPCW(CW) _MM_SET_ROUNDING_MODE(CW) #define CGAL_IA_GETFPCW(CW) CW = _MM_GET_ROUNDING_MODE() @@ -430,7 +497,8 @@ typedef unsigned int FPU_CW_t; #define CGAL_FE_UPWARD FP_RND_RP #define CGAL_FE_DOWNWARD FP_RND_RM -#elif defined ( _MSC_VER ) +#elif defined ( _MSC_VER ) && !defined ( _M_ARM64 ) +// x86/x64 MSVC only: _controlfp_s is not available on ARM64. #if ( _MSC_VER < 1400) #define CGAL_IA_SETFPCW(CW) _controlfp (CW, _MCW_RC ) #define CGAL_IA_GETFPCW(CW) CW = _controlfp (0, 0 ) & _MCW_RC diff --git a/Number_types/include/CGAL/Interval_nt.h b/Number_types/include/CGAL/Interval_nt.h index f70312178b75..341c042c4718 100644 --- a/Number_types/include/CGAL/Interval_nt.h +++ b/Number_types/include/CGAL/Interval_nt.h @@ -74,6 +74,8 @@ class Interval_nt #ifndef CGAL_NO_ASSERTIONS # ifdef CGAL_USE_SSE2 : val(_mm_setr_pd(-1, 0)) +# elif defined CGAL_USE_NEON + : val(vcombine_f64(vdup_n_f64(-1.0), vdup_n_f64(0.0))) # else : _inf(-1), _sup(0) # endif @@ -146,6 +148,9 @@ class Interval_nt // This constructor should really be private, like the simd() function, but // that would mean a lot of new friends, so they are only undocumented. explicit Interval_nt(__m128d v) : val(v) {} +#elif defined CGAL_USE_NEON + // NEON: float64x2_t holds {-inf, sup} just like __m128d. + explicit Interval_nt(float64x2_t v) : val(v) {} #endif // Unchecked version for Lazy_rep in Lazy.h. @@ -153,6 +158,9 @@ class Interval_nt Interval_nt(double i, double s, no_check_t) #ifdef CGAL_USE_SSE2 : val(_mm_setr_pd(-i, s)) +#elif defined CGAL_USE_NEON + // vsetq_lane_f64 builds {-i, s} in a float64x2_t register. + : val(vsetq_lane_f64(s, vdupq_n_f64(-i), 1)) #else : _inf(-i), _sup(s) #endif @@ -179,6 +187,9 @@ class Interval_nt { #ifdef CGAL_USE_SSE2 return IA (swap_m128d(val)); +#elif defined CGAL_USE_NEON + // Swapping {-inf, sup} gives {sup, -inf} = {-(-inf), -(sup)} = negated interval. + return IA (swap_f64x2(val)); #else return IA (-sup(), -inf()); #endif @@ -199,6 +210,11 @@ class Interval_nt #ifdef CGAL_USE_SSE2 // Faster to answer yes, but slower to answer no. return _mm_movemask_pd (_mm_cmpneq_pd (val, d.val)) == 0; +#elif defined CGAL_USE_NEON + // vceqq_f64 returns a uint64x2_t with all-ones lanes where equal. + // vandq_u64 combines both lanes; vgetq_lane_u64 checks the result. + uint64x2_t eq = vceqq_f64(val, d.val); + return (vgetq_lane_u64(eq, 0) & vgetq_lane_u64(eq, 1)) != 0; #else return inf() == d.inf() && sup() == d.sup(); #endif @@ -211,6 +227,17 @@ class Interval_nt __m128d y = _mm_xor_pd ((-d).val, m); // {-ds,di} __m128d c = _mm_cmplt_pd (val, y); // {i>ds,s {-ds, di} + // vcltq_f64(val, y) checks {-inf < -ds, sup < di} i.e. {inf>ds, sup sup() || d.sup() < inf()); #endif @@ -220,6 +247,9 @@ class Interval_nt { #ifdef CGAL_USE_SSE2 return -_mm_cvtsd_f64(val); +#elif defined CGAL_USE_NEON + // Lane 0 holds -inf; negate to get inf. + return -vgetq_lane_f64(val, 0); #else return -_inf; #endif @@ -232,12 +262,17 @@ class Interval_nt // - it is too opaque // - it is a less likely CSE candidate // return _mm_cvtsd_f64(_mm_unpackhi_pd(val, val)); +#elif defined CGAL_USE_NEON + // Lane 1 holds sup directly. + return vgetq_lane_f64(val, 1); #else return _sup; #endif } #ifdef CGAL_USE_SSE2 __m128d simd() const { return val; } +#elif defined CGAL_USE_NEON + float64x2_t simd() const { return val; } #endif std::pair pair() const @@ -272,6 +307,9 @@ class Interval_nt // is free to access. #ifdef CGAL_USE_SSE2 __m128d val; +#elif defined CGAL_USE_NEON + // float64x2_t stores {-inf, sup} in lanes {0, 1}, mirroring the SSE2 layout. + float64x2_t val; #else double _inf, _sup; #endif @@ -470,6 +508,12 @@ class Interval_nt __m128d bb = IA_opacify128_weak(b.simd()); __m128d r = _mm_add_pd(aa, bb); return Interval_nt(IA_opacify128(r)); +#elif defined CGAL_USE_NEON + // vaddq_f64 adds both lanes simultaneously: {-ai+-bi, as+bs}. + float64x2_t aa = IA_opacify_neon(a.simd()); + float64x2_t bb = IA_opacify_neon_weak(b.simd()); + float64x2_t r = vaddq_f64(aa, bb); + return Interval_nt(IA_opacify_neon(r)); #else return Interval_nt (-CGAL_IA_ADD(-a.inf(), -b.inf()), CGAL_IA_ADD(a.sup(), b.sup())); @@ -507,7 +551,7 @@ class Interval_nt Interval_nt operator- (const Interval_nt &a, const Interval_nt & b) { -#ifdef CGAL_USE_SSE2 +#if defined CGAL_USE_SSE2 || defined CGAL_USE_NEON return a+-b; #else Internal_protector P; @@ -516,7 +560,7 @@ class Interval_nt #endif } -#ifdef CGAL_USE_SSE2 +#if defined CGAL_USE_SSE2 || defined CGAL_USE_NEON friend Interval_nt operator- (double a, const Interval_nt & b) @@ -536,8 +580,38 @@ class Interval_nt else if(CGAL_CST_TRUE(b.is_point())) return a * b.inf(); #endif - Internal_protector P; -#ifdef CGAL_USE_SSE2 + Internal_protector P; +#ifdef CGAL_USE_NEON + // NEON multiply for interval arithmetic. + // compute all four cross-products and take the component-wise max. + // Layout: val = {-i, s} so lane0=-inf, lane1=sup. + float64x2_t aa = IA_opacify_neon_weak(a.simd()); // {-ai, as} + float64x2_t bb = b.simd(); // {-bi, bs} + float64x2_t neg_zero = vdupq_n_f64(-0.0); + float64x2_t neg_zero1 = vdupq_n_f64(-0.0); // same mask + float64x2_t ax = swap_f64x2(aa); // {as, -ai} + float64x2_t ap = vreinterpretq_f64_u64( + veorq_u64(vreinterpretq_u64_f64(ax), + vreinterpretq_u64_f64(neg_zero1))); // {-as, ai} + // bz = {bi, bs}: flip sign of lane0 of bb + float64x2_t bz_mask = vcombine_f64(vdup_n_f64(-0.0), vdup_n_f64(0.0)); + float64x2_t bz = vreinterpretq_f64_u64( + veorq_u64(vreinterpretq_u64_f64(bb), + vreinterpretq_u64_f64(bz_mask))); // {bi, bs} + bz = IA_opacify_neon(bz); + float64x2_t c = swap_f64x2(bz); // {bs, bi} + float64x2_t big = IA::largest().simd(); + float64x2_t x1 = vmulq_f64(aa, bz); // {-ai*bi, as*bs} + float64x2_t x2 = vmulq_f64(aa, c); // {-ai*bs, as*bi} + x2 = vminq_f64(x2, big); + float64x2_t x3 = vmulq_f64(ap, bz); // {-as*bi, ai*bs} + float64x2_t x4 = vmulq_f64(ap, c); // {-as*bs, ai*bi} + x4 = vminq_f64(x4, big); + float64x2_t y1 = vmaxq_f64(x1, x2); + float64x2_t y2 = vmaxq_f64(x3, x4); + float64x2_t r = vmaxq_f64(y1, y2); + return IA(IA_opacify_neon(r)); +#elif defined CGAL_USE_SSE2 # if !defined __SSE4_1__ && !defined __AVX__ // Brutal, compute all products in all directions. // The actual winner (by a hair) on recent hardware before removing NaNs. @@ -713,6 +787,13 @@ class Interval_nt // larger than necessary, but is likely faster to produce. r = _mm_min_pd(r,largest().simd()); return IA(IA_opacify128(r)); +#elif defined CGAL_USE_NEON + // Broadcast scalar a into both lanes, multiply, clamp infinities. + float64x2_t bb = IA_opacify_neon_weak(b.simd()); + float64x2_t aa = vdupq_n_f64(IA_opacify(a)); + float64x2_t r = vmulq_f64(aa, bb); + r = vminq_f64(r, largest().simd()); + return IA(IA_opacify_neon(r)); #else else if (!(a > 0)) return 0.; // We could test this before the SSE block and remove the minpd line. return IA(-CGAL_IA_MUL(a, -b.inf()), CGAL_IA_MUL(a, b.sup())); @@ -737,8 +818,36 @@ class Interval_nt else if(CGAL_CST_TRUE(b.is_point())) return a / b.inf(); #endif - Internal_protector P; -#if defined CGAL_USE_SSE2 && (defined __SSE4_1__ || defined __AVX__) + Internal_protector P; +#if defined CGAL_USE_NEON + // Division for interval arithmetic on NEON. + // Check whether 0 is in b: if both -bi<=0 and bs>=0 then bi<=0 && bs>=0. + // val layout {-bi, bs}: lane0=-bi, lane1=bs. + // 0 in b <=> -bi <= 0 AND bs >= 0 <=> lane0 <= 0 AND lane1 >= 0. + float64x2_t zero = vdupq_n_f64(0.0); + uint64x2_t ge0 = vcgeq_f64(b.simd(), zero); // {-bi>=0, bs>=0} + if (vgetq_lane_u64(ge0, 0) && vgetq_lane_u64(ge0, 1)) + return largest(); // 0 in b + // b is strictly positive or strictly negative; use multiply-by-reciprocal + // approach: compute {1/-bi, 1/bs} then multiply. + // For b>0: {-bi,bs} both positive -> reciprocal is safe. + // For b<0: negate b first (swap lanes), divide, negate result. + float64x2_t bb = b.simd(); + float64x2_t aa = a.simd(); + // Determine sign of b from lane1 (bs): if bs < 0 then b < 0. + if (vgetq_lane_f64(bb, 1) < 0.0) { + // b < 0: negate both a and b, then divide as if b > 0. + aa = swap_f64x2(aa); // negate a + bb = swap_f64x2(bb); // negate b + } + // Now b > 0: {-bi, bs} with bi < 0 < bs. + // Compute reciprocal of b lanes and multiply. + float64x2_t ones = vdupq_n_f64(1.0); + float64x2_t inv_b = vdivq_f64(ones, IA_opacify_neon(bb)); + float64x2_t r = vdivq_f64(IA_opacify_neon_weak(aa), + IA_opacify_neon(bb)); + return IA(IA_opacify_neon(r)); +#elif defined CGAL_USE_SSE2 && (defined __SSE4_1__ || defined __AVX__) //// not a tight bound, but easy: // return CGAL::inverse(b)*a; # if 1 @@ -859,6 +968,43 @@ class Interval_nt __m128d r = _mm_div_pd(aa, bb); return Interval_nt(IA_opacify128(r)); } +#elif defined CGAL_USE_NEON + friend + Interval_nt + operator/ (double a, const Interval_nt & b) + { + // Check 0 in b: lane0=-bi>=0 AND lane1=bs>=0. + float64x2_t zero = vdupq_n_f64(0.0); + uint64x2_t ge0 = vcgeq_f64(b.simd(), zero); + if (vgetq_lane_u64(ge0, 0) && vgetq_lane_u64(ge0, 1)) + return largest(); + float64x2_t aa, xx; + if (a > 0) { + aa = vdupq_n_f64(-a); + xx = (-b).simd(); + } else if (a < 0) { + aa = vdupq_n_f64(a); + xx = b.simd(); + } else return 0.; + Internal_protector P; + float64x2_t r = vdivq_f64(IA_opacify_neon_weak(aa), IA_opacify_neon(xx)); + return Interval_nt(IA_opacify_neon(r)); + } + + friend + Interval_nt + operator/ (Interval_nt a, double b) + { + if (b < 0) { a = -a; b = -b; } + else if (b == 0) return largest(); + // Now b > 0 + Internal_protector P; + b = IA_opacify(b); + float64x2_t bb = vdupq_n_f64(b); + float64x2_t aa = IA_opacify_neon(a.simd()); + float64x2_t r = vdivq_f64(aa, bb); + return Interval_nt(IA_opacify_neon(r)); + } #endif }; @@ -903,6 +1049,15 @@ magnitude (const Interval_nt & d) __m128d x = _mm_and_pd (d.simd(), m); // { abs(inf), abs(sup) } __m128d y = _mm_unpackhi_pd (x, x); return _mm_cvtsd_f64 (_mm_max_sd (x, y)); +#elif defined CGAL_USE_NEON + // Mask off sign bits to get absolute values of both lanes. + // val = {-inf, sup}; abs gives {|inf|, |sup|}; return the max. + const uint64x2_t abs_mask = vdupq_n_u64(0x7fffffffffffffffULL); + float64x2_t x = vreinterpretq_f64_u64( + vandq_u64(vreinterpretq_u64_f64(d.simd()), abs_mask)); + float64_t lo = vgetq_lane_f64(x, 0); + float64_t hi = vgetq_lane_f64(x, 1); + return lo > hi ? lo : hi; #else return (std::max)(CGAL::abs(d.inf()), CGAL::abs(d.sup())); #endif @@ -978,6 +1133,17 @@ struct Min > // Use _mm_max_sd instead? __m128d y = _mm_max_pd (d.simd(), e.simd()); return Interval_nt (_mm_move_sd (x, y)); +#elif defined CGAL_USE_NEON + // min({-di,ds},{-ei,es}) = {-max(di,ei), min(ds,es)} + // lane0: min(-di,-ei) = -max(di,ei) -> vminq_f64 on lane0 + // lane1: min(ds,es) -> vminq_f64 on lane1 + // But we want lane0 from vmax (most-negative = largest -inf) + // and lane1 from vmin (smallest sup). + float64x2_t x = vminq_f64(d.simd(), e.simd()); // {min(-di,-ei), min(ds,es)} + float64x2_t y = vmaxq_f64(d.simd(), e.simd()); // {max(-di,-ei), max(ds,es)} + // Result: lane0 from y (largest -inf = most negative lower bound), + // lane1 from x (smallest sup). + return Interval_nt(vcombine_f64(vget_low_f64(y), vget_high_f64(x))); #else return Interval_nt( -(std::max)(-d.inf(), -e.inf()), @@ -1000,6 +1166,15 @@ struct Max > __m128d x = _mm_min_pd (d.simd(), e.simd()); __m128d y = _mm_max_pd (d.simd(), e.simd()); return Interval_nt (_mm_move_sd (y, x)); +#elif defined CGAL_USE_NEON + // max({-di,ds},{-ei,es}) = {-min(di,ei), max(ds,es)} + // lane0: min(-di,-ei) = -max(di,ei) -> vminq_f64 gives smallest -inf + // lane1: max(ds,es) -> vmaxq_f64 + float64x2_t x = vminq_f64(d.simd(), e.simd()); // {min(-di,-ei), min(ds,es)} + float64x2_t y = vmaxq_f64(d.simd(), e.simd()); // {max(-di,-ei), max(ds,es)} + // Result: lane0 from x (smallest -inf = least negative lower bound), + // lane1 from y (largest sup). + return Interval_nt(vcombine_f64(vget_low_f64(x), vget_high_f64(y))); #else return Interval_nt( -(std::min)(-d.inf(), -e.inf()), @@ -1132,7 +1307,16 @@ namespace INTERN_INTERVAL_NT { // sqrt([+a,+b]) => [sqrt(+a);sqrt(+b)] // sqrt([-a,+b]) => [0;sqrt(+b)] => assumes roundoff error. // sqrt([-a,-b]) => [0;sqrt(-b)] => assumes user bug (unspecified result). -#ifdef __AVX512F__ +#ifdef CGAL_USE_NEON + // AArch64 vsqrtq_f64 computes sqrt in round-to-nearest mode. + // We need the lower bound rounded DOWN and the upper bound rounded UP. + // Since we cannot change rounding mode per-lane, we use the same + // nextafter-based approach as CGAL_ALWAYS_ROUND_TO_NEAREST. + double i = 0.0; + if (d.inf() > 0.0) + i = nextafter(std::sqrt(d.inf()), 0.0); + return Interval_nt(i, IA_sqrt_up(d.sup())); +#elif defined __AVX512F__ double i = 0; if(d.inf() > 0){ __m128d x = d.simd(); @@ -1171,6 +1355,16 @@ namespace INTERN_INTERVAL_NT { __m128d b = _mm_xor_pd(a, _mm_setr_pd(-0., 0.)); // {i,s} __m128d r = _mm_mul_pd(a, b); // {-i*i,s*s} return Interval_nt(IA_opacify128(r)); +#elif defined CGAL_USE_NEON + // abs(d).simd() = {-i, s} with 0 <= i <= s. + // Flip sign of lane0 to get {i, s}, then multiply: {-i*i, s*s}. + float64x2_t a = IA_opacify_neon(CGAL::abs(d).simd()); // {-i, s} + float64x2_t sign_flip = vcombine_f64(vdup_n_f64(-0.0), vdup_n_f64(0.0)); + float64x2_t b = vreinterpretq_f64_u64( + veorq_u64(vreinterpretq_u64_f64(a), + vreinterpretq_u64_f64(sign_flip))); // {i, s} + float64x2_t r = vmulq_f64(a, b); // {-i*i, s*s} + return Interval_nt(IA_opacify_neon(r)); #else if (d.inf()>=0.0) return Interval_nt(-CGAL_IA_MUL(-d.inf(), d.inf()), @@ -1197,6 +1391,21 @@ namespace INTERN_INTERVAL_NT { __m128d z = _mm_set1_pd(-0.); // +0. would be valid, but I'd rather end up with interval [+0, sup] __m128d r = _mm_min_sd(t, z); return Interval_nt (r); +#elif defined CGAL_USE_NEON + // a = {-inf, sup}, -a = swap = {sup, -inf} = {-(-inf), -(sup)} + // min lane-wise gives the more-negative lower bound (larger -inf). + // max lane-wise gives the larger upper bound. + // Then combine: lane0 from min (clamped to <=0), lane1 from max. + float64x2_t a = d.simd(); + float64x2_t b = (-d).simd(); + float64x2_t x = vminq_f64(a, b); + float64x2_t y = vmaxq_f64(a, b); + // Combine: take lane0 from x (lower bound side), lane1 from y (upper bound). + float64x2_t t = vcombine_f64(vget_low_f64(x), vget_high_f64(y)); + // Clamp lane0 to <= 0 (same as _mm_min_sd with -0.). + float64x2_t neg_zero_lo = vcombine_f64(vdup_n_f64(-0.0), vdup_n_f64(0.0)); + float64x2_t r = vminq_f64(t, neg_zero_lo); + return Interval_nt(r); #else if (d.inf() >= 0.0) return d; if (d.sup() <= 0.0) return -d; @@ -1548,6 +1757,10 @@ class Interval_traits< Interval_nt > Interval operator()( const Interval& a, const Interval& b ) const { #ifdef CGAL_USE_SSE2 return Interval(_mm_max_pd(a.simd(), b.simd())); +#elif defined CGAL_USE_NEON + // vmaxq_f64 on {-ai,as} and {-bi,bs} gives {max(-ai,-bi), max(as,bs)} + // = {-min(ai,bi), max(as,bs)} which is exactly the hull. + return Interval(vmaxq_f64(a.simd(), b.simd())); #else BOOST_USING_STD_MAX(); BOOST_USING_STD_MIN(); diff --git a/Number_types/include/CGAL/double.h b/Number_types/include/CGAL/double.h index b673205646c6..deaaa41dd1eb 100644 --- a/Number_types/include/CGAL/double.h +++ b/Number_types/include/CGAL/double.h @@ -30,6 +30,11 @@ #include #endif +// NEON fabs: enabled when CGAL_USE_NEON_FABS is defined. +#ifdef CGAL_USE_NEON_FABS +# include +#endif + #ifdef _MSC_VER #include #endif @@ -132,7 +137,20 @@ inline double sse2fabs(double a) temp = _mm_and_pd(temp, absMask.m); return _mm_cvtsd_f64 (temp); } +#endif +// NEON absolute value for double. +#ifdef CGAL_USE_NEON_FABS +inline double neon_fabs(double a) +{ +#if defined(_MSC_VER) && defined(_M_ARM64) + // MSVC ARM64 does not provide vabsd_f64(). + return vget_lane_f64(vabs_f64(vdup_n_f64(a)), 0); +#else + // GCC and Clang map this to a single FABS instruction. + return vabsd_f64(a); +#endif +} #endif template <> class Real_embeddable_traits< double > @@ -141,13 +159,15 @@ template <> class Real_embeddable_traits< double > // GCC is faster with std::fabs(). -#if defined(__GNUG__) || defined(CGAL_MSVC_USE_STD_FABS) || defined(CGAL_USE_SSE2_FABS) +#if defined(__GNUG__) || defined(CGAL_MSVC_USE_STD_FABS) || defined(CGAL_USE_SSE2_FABS) || defined(CGAL_USE_NEON_FABS) class Abs : public CGAL::cpp98::unary_function< Type, Type > { public: Type operator()( const Type& x ) const { #ifdef CGAL_USE_SSE2_FABS return sse2fabs(x); +#elif defined CGAL_USE_NEON_FABS + return neon_fabs(x); #else return std::fabs( x ); #endif diff --git a/Number_types/include/CGAL/sse2.h b/Number_types/include/CGAL/sse2.h index 27a56376dd62..e41f7d1e1493 100644 --- a/Number_types/include/CGAL/sse2.h +++ b/Number_types/include/CGAL/sse2.h @@ -17,12 +17,18 @@ #ifndef CGAL_SSE2_H #define CGAL_SSE2_H -#include #if defined ( _MSC_VER ) #define CGAL_ALIGN_16 __declspec(align(16)) -#elif defined( __GNUC__ ) +#elif defined( __GNUC__ ) || defined(__clang__) #define CGAL_ALIGN_16 __attribute__((aligned(16))) #endif +// Architecture-specific SIMD headers +#if defined(__aarch64__) || defined(_M_ARM64) +# include +#else +# include +#endif + #endif // CGAL_SSE2_H diff --git a/Number_types/include/CGAL/utils_classes.h b/Number_types/include/CGAL/utils_classes.h index f7052a8c59a0..22cc4c18d982 100644 --- a/Number_types/include/CGAL/utils_classes.h +++ b/Number_types/include/CGAL/utils_classes.h @@ -21,6 +21,12 @@ #include #endif +// NEON min/max for double: enabled by CGAL_USE_NEON_MAX. +// vmaxq_f64 / vminq_f64 operate on two doubles in a single 128-bit register. +#ifdef CGAL_USE_NEON_MAX +# include +#endif + namespace CGAL { template < class A, class B = A > @@ -194,14 +200,95 @@ inline void sse2minmax(double& a, double b, double& c) #endif // CGAL_USE_SSE2_MAX +#ifdef CGAL_USE_NEON_MAX + +inline double neon_max(double a, double b) +{ + float64x2_t A = vdupq_n_f64(a); + float64x2_t B = vdupq_n_f64(b); + return vgetq_lane_f64(vmaxq_f64(A, B), 0); +} + +inline double neon_max(double a, double b, double c) +{ + float64x2_t AB = vmaxq_f64(vdupq_n_f64(a), vdupq_n_f64(b)); + float64x2_t C = vdupq_n_f64(c); + return vgetq_lane_f64(vmaxq_f64(AB, C), 0); +} + +inline double neon_max(double a, double b, double c, double d) +{ + float64x2_t AB = vmaxq_f64(vdupq_n_f64(a), vdupq_n_f64(b)); + float64x2_t CD = vmaxq_f64(vdupq_n_f64(c), vdupq_n_f64(d)); + return vgetq_lane_f64(vmaxq_f64(AB, CD), 0); +} + +inline double neon_min(double a, double b) +{ + float64x2_t A = vdupq_n_f64(a); + float64x2_t B = vdupq_n_f64(b); + return vgetq_lane_f64(vminq_f64(A, B), 0); +} + +inline double neon_min(double a, double b, double c) +{ + float64x2_t AB = vminq_f64(vdupq_n_f64(a), vdupq_n_f64(b)); + float64x2_t C = vdupq_n_f64(c); + return vgetq_lane_f64(vminq_f64(AB, C), 0); +} + +inline double neon_min(double a, double b, double c, double d) +{ + float64x2_t AB = vminq_f64(vdupq_n_f64(a), vdupq_n_f64(b)); + float64x2_t CD = vminq_f64(vdupq_n_f64(c), vdupq_n_f64(d)); + return vgetq_lane_f64(vminq_f64(AB, CD), 0); +} + +// sets a = min(a,b,c), c = max(a,b,c). +// b may hold any value on exit (same contract as sse2minmax). +inline void neon_minmax(double& a, double b, double& c) +{ + float64x2_t A = vdupq_n_f64(a); + float64x2_t B = vdupq_n_f64(b); + float64x2_t C = vdupq_n_f64(c); + + float64x2_t AB_min = vminq_f64(A, B); + float64x2_t AB_max = vmaxq_f64(A, B); + + // min(a,b,c) + float64x2_t min_abc = vminq_f64(AB_min, C); + // max(a,b,c) + float64x2_t max_abc = vmaxq_f64(AB_max, C); + + a = vgetq_lane_f64(min_abc, 0); + c = vgetq_lane_f64(max_abc, 0); +} + +#endif // CGAL_USE_NEON_MAX + +// simd_minmax: unified name for call sites that work on both SSE2 and NEON. +// Defined only when one of the two SIMD back-ends is active. +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX +inline void simd_minmax(double& a, double b, double& c) +{ +#ifdef CGAL_USE_SSE2_MAX + sse2minmax(a, b, c); +#else + neon_minmax(a, b, c); +#endif +} +#endif // CGAL_USE_SSE2_MAX || CGAL_USE_NEON_MAX + template <> struct Max :public CGAL::cpp98::binary_function< double, double, double > { - Max() {} +Max() {} - double operator()( const double& x, const double& y) const +double operator()( const double& x, const double& y) const { #ifdef CGAL_USE_SSE2_MAX return sse2max(x,y); +#elif defined CGAL_USE_NEON_MAX + return neon_max(x,y); #else return (std::max)( x, y); #endif @@ -211,6 +298,8 @@ struct Max :public CGAL::cpp98::binary_function< double, double, double { #ifdef CGAL_USE_SSE2_MAX return sse2max(x,y,z); +#elif defined CGAL_USE_NEON_MAX + return neon_max(x,y,z); #else return (std::max)((std::max)( x, y), z); #endif @@ -220,6 +309,8 @@ struct Max :public CGAL::cpp98::binary_function< double, double, double { #ifdef CGAL_USE_SSE2_MAX return sse2max(w,x,y,z); +#elif defined CGAL_USE_NEON_MAX + return neon_max(w,x,y,z); #else return (std::max)((std::max)( x, y), (std::max)(w,z)); #endif @@ -234,6 +325,8 @@ struct Min :public CGAL::cpp98::binary_function< double, double, double { #ifdef CGAL_USE_SSE2_MAX return sse2min(x,y); +#elif defined CGAL_USE_NEON_MAX + return neon_min(x,y); #else return (std::min)( x, y); #endif @@ -243,6 +336,8 @@ struct Min :public CGAL::cpp98::binary_function< double, double, double { #ifdef CGAL_USE_SSE2_MAX return sse2min(x,y,z); +#elif defined CGAL_USE_NEON_MAX + return neon_min(x,y,z); #else return (std::min)((std::min)( x, y), z); #endif @@ -252,6 +347,8 @@ struct Min :public CGAL::cpp98::binary_function< double, double, double { #ifdef CGAL_USE_SSE2_MAX return sse2min(w,x,y,z); +#elif defined CGAL_USE_NEON_MAX + return neon_min(w,x,y,z); #else return (std::min)((std::min)( x, y), (std::min)(w,z)); #endif diff --git a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Static_filters/Periodic_3_orientation_3.h b/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Static_filters/Periodic_3_orientation_3.h index f858558c45f7..652fa5abbd4d 100644 --- a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Static_filters/Periodic_3_orientation_3.h +++ b/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Static_filters/Periodic_3_orientation_3.h @@ -100,7 +100,7 @@ class Periodic_3_orientation_3 double aprz = CGAL::abs(prz); double apsz = CGAL::abs(psz); -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX CGAL::Max mmax; maxx = mmax(maxx, aprx, apsx); @@ -121,14 +121,14 @@ class Periodic_3_orientation_3 double eps = 5.1107127829973299e-15 * maxx * maxy * maxz; -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX /* CGAL::Min mmin; double tmp = mmin(maxx, maxy, maxz); maxz = mmax(maxx, maxy, maxz); maxx = tmp; */ - sse2minmax(maxx,maxy,maxz); + simd_minmax(maxx,maxy,maxz); // maxy can contain ANY element #else // Sort maxx < maxy < maxz. @@ -218,7 +218,7 @@ class Periodic_3_orientation_3 double apsy = CGAL::abs(psy); double apsz = CGAL::abs(psz); -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX CGAL::Max mmax; maxx = mmax(maxx, aqtx, artx, astx); maxy = mmax(maxy, aqty, arty, asty); @@ -239,14 +239,14 @@ class Periodic_3_orientation_3 double eps = 4.111024169857068197e-15 * maxx * maxy * maxz; -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX /* CGAL::Min mmin; double tmp = mmin(maxx, maxy, maxz); maxz = mmax(maxx, maxy, maxz); maxx = tmp; */ - sse2minmax(maxx,maxy,maxz); + simd_minmax(maxx,maxy,maxz); // maxy can contain ANY element #else // Sort maxx < maxy < maxz. diff --git a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Static_filters/Periodic_3_side_of_oriented_sphere_3.h b/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Static_filters/Periodic_3_side_of_oriented_sphere_3.h index 9c0eae4b8de8..7f9e4850729b 100644 --- a/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Static_filters/Periodic_3_side_of_oriented_sphere_3.h +++ b/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Static_filters/Periodic_3_side_of_oriented_sphere_3.h @@ -113,7 +113,7 @@ class Periodic_3_side_of_oriented_sphere_3 double asty = CGAL::abs(sty); double astz = CGAL::abs(stz); -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX CGAL::Max mmax; maxx = mmax(maxx, aqtx, artx, astx); maxy = mmax(maxy, aqty, arty, asty); @@ -134,14 +134,14 @@ class Periodic_3_side_of_oriented_sphere_3 double eps = 1.2466136531027298e-13 * maxx * maxy * maxz; -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX /* CGAL::Min mmin; double tmp = mmin(maxx, maxy, maxz); maxz = mmax(maxx, maxy, maxz); maxx = tmp; */ - sse2minmax(maxx,maxy,maxz); + simd_minmax(maxx,maxy,maxz); // maxy can contain ANY element #else @@ -255,7 +255,7 @@ class Periodic_3_side_of_oriented_sphere_3 double asty = CGAL::abs(sty); double astz = CGAL::abs(stz); -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX CGAL::Max mmax; maxx = mmax(maxx, aqtx, artx, astx); maxy = mmax(maxy, aqty, arty, asty); @@ -276,14 +276,14 @@ class Periodic_3_side_of_oriented_sphere_3 double eps = 1.0466759304746772485e-13 * maxx * maxy * maxz; -#ifdef CGAL_USE_SSE2_MAX +#if defined CGAL_USE_SSE2_MAX || defined CGAL_USE_NEON_MAX /* CGAL::Min mmin; double tmp = mmin(maxx, maxy, maxz); maxz = mmax(maxx, maxy, maxz); maxx = tmp; */ - sse2minmax(maxx,maxy,maxz); + simd_minmax(maxx,maxy,maxz); // maxy can contain ANY element #else // Sort maxx < maxy < maxz.