diff --git a/.github/workflows/setup.yml b/.github/workflows/setup.yml new file mode 100755 index 0000000..c3e06e6 --- /dev/null +++ b/.github/workflows/setup.yml @@ -0,0 +1,39 @@ +name: predicates test + +on: + workflow_dispatch: + pull_request: + push: + +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest] + + steps: + - uses: actions/checkout@v3 + + - name: Print environment + run: | + echo github.event.action: ${{ github.event.action }} + echo github.event_name: ${{ github.event_name }} + + - name: Install dependencies on unx + if: startsWith(matrix.os, 'ubuntu') + run: | + sudo apt-get update + g++ --version + + - name: Build predicates + run: | + cd ${{github.workspace}} + ls + g++ -std=c++17 -pedantic -Wall -O3 -flto -DNDEBUG \ + example.cpp -oexample + + - name: Check predicates + run: | + ./example diff --git a/README.md b/README.md index ad45de8..0dbafe1 100644 --- a/README.md +++ b/README.md @@ -19,6 +19,10 @@ bisect3d: orientation of point wrt. half-space in E^3. bisect3w: orientation of point wrt. half-space in E^3 (weighted). inball3d: point-in-circumball (Delaunay-Voronoi tessellations) in E^3. inball3w: point-in-ortho-ball (Regular-Laguerre tessellations) in E^3. + +orient4d: orientation of 5 points in E^4, or a point wrt. a cell. +bisect4d: orientation of point wrt. half-space in E^4. +bisect4w: orientation of point wrt. half-space in E^4 (weighted). ```` A simplified two-stage variation on Shewchuk's original arithmetic is employed, adopting standard (fast!) floating-point approximations when results are unambiguous and falling back onto (slower) arbitrary precision evaluations as necessary to guarantee "sign-correctness". Semi-static filters are used to toggle between floating-point and arbitrary precision kernels. An optional third (intermediate) stage based on interval arithmetic is also available. diff --git a/basebase.hpp b/basebase.hpp index ecd4e1c..f99e6df 100644 --- a/basebase.hpp +++ b/basebase.hpp @@ -31,11 +31,11 @@ * ------------------------------------------------------------ * - * Last updated: 02 March, 2020 + * Last updated: 11 May., 2024 * - * Copyright 2013-2020 + * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * ------------------------------------------------------------ diff --git a/example.cpp b/example.cpp index df4c37d..52ff2f7 100644 --- a/example.cpp +++ b/example.cpp @@ -176,6 +176,62 @@ int main () { std::cout << "inball3w: " << _rr; std::cout << std::endl; +/*-------------------------------- test prediactes in E^4 */ + + double _P1[5] = { // (d+1) coord. is weight + +0.0, +0.0, +0.0, +0.0, +0.0 + } ; + double _P2[5] = { + +1.0, +0.0, +0.0, +0.0, +0.1 + } ; + double _P3[5] = { + +1.0, +1.0, +0.0, +0.0, +0.2 + } ; + double _P4[5] = { + +1.0, +1.0, +1.0, +0.0, +0.3 + } ; +// double _P5[5] = { +// +1.0, +1.0, +1.0, +1.0, +0.4 +// } ; + + double _Q4[5] = { + +0.5, +0.5, +0.5, +0.5, +0.0 + } ; + + // Test the orienation of the point QQ wrt. the cell + // P1, P2, P3, P4 in E^4. + + _rr = geompred::orient4d ( + _P1, _P2, _P3, _P4, _Q4 + ) ; + + std::cout << "orient4d: " << _rr; + std::cout << std::endl; + + // Test the orienation of the point QQ wrt. the half- + // space of P1, P2. + + // This is the unweighted case in E^4. + + _rr = geompred::bisect4d ( + _P1, _P2, _Q4 + ) ; + + std::cout << "bisect4d: " << _rr; + std::cout << std::endl; + + // Test the orienation of the point QQ wrt. the half- + // space of P1, P2. + + // This is the "weighted" case in E^4. + + _rr = geompred::bisect4w ( + _P1, _P2, _Q4 + ) ; + + std::cout << "bisect4w: " << _rr; + std::cout << std::endl; + return 0 ; } diff --git a/expansion/dd_float.hpp b/expansion/dd_float.hpp index 7b13176..c94bd39 100644 --- a/expansion/dd_float.hpp +++ b/expansion/dd_float.hpp @@ -37,20 +37,24 @@ * how they can obtain it for free, then you are not * required to make any arrangement with me.) * - * Disclaimer: Neither I nor: Columbia University, The - * Massachusetts Institute of Technology, The - * University of Sydney, nor The National Aeronautics - * and Space Administration warrant this code in any - * way whatsoever. This code is provided "as-is" to be - * used at your own risk. + * Disclaimer: Neither I nor THE CONTRIBUTORS warrant + * this code in any way whatsoever. This code is + * provided "as-is" to be used at your own risk. + * + * THE CONTRIBUTORS include: + * (a) The University of Sydney + * (b) The Massachusetts Institute of Technology + * (c) Columbia University + * (d) The National Aeronautics & Space Administration + * (e) Los Alamos National Laboratory * -------------------------------------------------------- * - * Last updated: 16 April, 2020 + * Last updated: 16 Apr., 2020 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * -------------------------------------------------------- diff --git a/expansion/ia_float.hpp b/expansion/ia_float.hpp index 14e9d0c..1d2cc80 100644 --- a/expansion/ia_float.hpp +++ b/expansion/ia_float.hpp @@ -22,20 +22,24 @@ * how they can obtain it for free, then you are not * required to make any arrangement with me.) * - * Disclaimer: Neither I nor: Columbia University, The - * Massachusetts Institute of Technology, The - * University of Sydney, nor The National Aeronautics - * and Space Administration warrant this code in any - * way whatsoever. This code is provided "as-is" to be - * used at your own risk. + * Disclaimer: Neither I nor THE CONTRIBUTORS warrant + * this code in any way whatsoever. This code is + * provided "as-is" to be used at your own risk. + * + * THE CONTRIBUTORS include: + * (a) The University of Sydney + * (b) The Massachusetts Institute of Technology + * (c) Columbia University + * (d) The National Aeronautics & Space Administration + * (e) Los Alamos National Laboratory * -------------------------------------------------------- * - * Last updated: 10 April, 2020 + * Last updated: 10 Apr., 2020 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * -------------------------------------------------------- diff --git a/expansion/mp_basic.hpp b/expansion/mp_basic.hpp index 3076ca1..66262e5 100644 --- a/expansion/mp_basic.hpp +++ b/expansion/mp_basic.hpp @@ -50,20 +50,24 @@ * how they can obtain it for free, then you are not * required to make any arrangement with me.) * - * Disclaimer: Neither I nor: Columbia University, The - * Massachusetts Institute of Technology, The - * University of Sydney, nor The National Aeronautics - * and Space Administration warrant this code in any - * way whatsoever. This code is provided "as-is" to be - * used at your own risk. + * Disclaimer: Neither I nor THE CONTRIBUTORS warrant + * this code in any way whatsoever. This code is + * provided "as-is" to be used at your own risk. + * + * THE CONTRIBUTORS include: + * (a) The University of Sydney + * (b) The Massachusetts Institute of Technology + * (c) Columbia University + * (d) The National Aeronautics & Space Administration + * (e) Los Alamos National Laboratory * -------------------------------------------------------- * - * Last updated: 16 April, 2020 + * Last updated: 16 Apr., 2020 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * -------------------------------------------------------- diff --git a/expansion/mp_float.hpp b/expansion/mp_float.hpp index ca7a7ef..a38b452 100644 --- a/expansion/mp_float.hpp +++ b/expansion/mp_float.hpp @@ -46,20 +46,24 @@ * how they can obtain it for free, then you are not * required to make any arrangement with me.) * - * Disclaimer: Neither I nor: Columbia University, The - * Massachusetts Institute of Technology, The - * University of Sydney, nor The National Aeronautics - * and Space Administration warrant this code in any - * way whatsoever. This code is provided "as-is" to be - * used at your own risk. + * Disclaimer: Neither I nor THE CONTRIBUTORS warrant + * this code in any way whatsoever. This code is + * provided "as-is" to be used at your own risk. + * + * THE CONTRIBUTORS include: + * (a) The University of Sydney + * (b) The Massachusetts Institute of Technology + * (c) Columbia University + * (d) The National Aeronautics & Space Administration + * (e) Los Alamos National Laboratory * -------------------------------------------------------- * - * Last updated: 07 April, 2020 + * Last updated: 05 Jun., 2022 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * -------------------------------------------------------- @@ -95,7 +99,7 @@ indx_type static constexpr _size = N ; - real_type _xdat [ N ] ; + real_type _xdat [ N ] = { 0. } ; indx_type _xlen = 0 ; public : @@ -930,7 +934,8 @@ template < size_t AX, size_t BX, size_t AY, - size_t BY, size_t NP + size_t BY, + size_t NP > __inline_call void expansion_dot ( expansion const& _xa, @@ -976,6 +981,38 @@ expansion_add(_xp, _yp, _zp, _dp); } + template < + size_t AX, size_t BX, size_t AY, + size_t BY, size_t AZ, size_t BZ, + size_t AQ, size_t BQ, + size_t NP + > + __inline_call void expansion_dot ( + expansion const& _xa, + expansion const& _xb, + expansion const& _ya, + expansion const& _yb, + expansion const& _za, + expansion const& _zb, + expansion const& _qa, + expansion const& _qb, + expansion & _dp + ) // 4-dim dotproduct + { + expansion _xp ; + expansion_mul(_xa, _xb, _xp); + + expansion _yp ; + expansion_mul(_ya, _yb, _yp); + + expansion _zp ; + expansion_mul(_za, _zb, _zp); + + expansion _qp ; + expansion_mul(_qa, _qb, _qp); + + expansion_add(_xp, _yp, _zp, _qp, _dp) ; + } # undef REAL_TYPE # undef INDX_TYPE diff --git a/expansion/mp_utils.hpp b/expansion/mp_utils.hpp index d35eed8..4f4e2aa 100644 --- a/expansion/mp_utils.hpp +++ b/expansion/mp_utils.hpp @@ -22,20 +22,24 @@ * how they can obtain it for free, then you are not * required to make any arrangement with me.) * - * Disclaimer: Neither I nor: Columbia University, The - * Massachusetts Institute of Technology, The - * University of Sydney, nor The National Aeronautics - * and Space Administration warrant this code in any - * way whatsoever. This code is provided "as-is" to be - * used at your own risk. + * Disclaimer: Neither I nor THE CONTRIBUTORS warrant + * this code in any way whatsoever. This code is + * provided "as-is" to be used at your own risk. + * + * THE CONTRIBUTORS include: + * (a) The University of Sydney + * (b) The Massachusetts Institute of Technology + * (c) Columbia University + * (d) The National Aeronautics & Space Administration + * (e) Los Alamos National Laboratory * -------------------------------------------------------- * - * Last updated: 03 March, 2020 + * Last updated: 11 May., 2024 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * -------------------------------------------------------- @@ -284,8 +288,8 @@ } else { - expansion_sub(_det2p, _det1p, _sum_1); - expansion_sub(_det4p, _det3p, _sum_2); + expansion_sub(_det1p, _det2p, _sum_1); + expansion_sub(_det3p, _det4p, _sum_2); } expansion_add(_sum_1, _sum_2, _final); @@ -436,6 +440,163 @@ expansion_add(_sum_3, _sum_2, _final); } + /* + -------------------------------------------------------- + * + * Compute an exact 6 x 6 determinant. + * + * | a1 a2 a3 a4 a5 v1 | + * | b1 b2 b3 b4 b5 v2 | + * | c1 c2 c3 c4 c5 v3 | + * | d1 d2 d3 d4 d5 v4 | + * | e1 e2 e3 e4 e5 v5 | + * | f1 f2 f3 f4 f5 v6 | + * + * as the product of 5 x 5 minors about a pivot column + * P, shown here for P = 6. The entry V1 is associated + * with the minor + * + * | b1 b2 b3 b4 b5 | + * | c1 c2 c3 c4 c5 | + * | d1 d2 d3 d4 d5 | = D1 + * | e1 e2 e3 e4 e5 | + * | f1 f2 f3 f4 f5 | + * + * and so on for (V2,D2), (V3,D3) etc. + * + -------------------------------------------------------- + */ + + template < + size_t NA, size_t NB, size_t NC, + size_t ND, size_t NE, size_t NF, + size_t NG, size_t NH, size_t NI, + size_t NJ, size_t NK, size_t NL, + size_t NM + > + __inline_call void compute_det_6x6 ( + expansion const& _det1p , + expansion const& _val1p , + expansion const& _det2p , + expansion const& _val2p , + expansion const& _det3p , + expansion const& _val3p , + expansion const& _det4p , + expansion const& _val4p , + expansion const& _det5p , + expansion const& _val5p , + expansion const& _det6p , + expansion const& _val6p , + expansion & _final , + INDX_TYPE _pivot + ) + { + /*---------------------------------- products Vi * Di */ + INDX_TYPE + constexpr N1 = mul_alloc (NA, NB) ; + expansion _mul1p; + expansion_mul(_det1p, _val1p, _mul1p); + + INDX_TYPE + constexpr N2 = mul_alloc (NC, ND) ; + expansion _mul2p; + expansion_mul(_det2p, _val2p, _mul2p); + + INDX_TYPE + constexpr N3 = mul_alloc (NE, NF) ; + expansion _mul3p; + expansion_mul(_det3p, _val3p, _mul3p); + + INDX_TYPE + constexpr N4 = mul_alloc (NG, NH) ; + expansion _mul4p; + expansion_mul(_det4p, _val4p, _mul4p); + + INDX_TYPE + constexpr N5 = mul_alloc (NI, NJ) ; + expansion _mul5p; + expansion_mul(_det5p, _val5p, _mul5p); + + INDX_TYPE + constexpr N6 = mul_alloc (NK, NL) ; + expansion _mul6p; + expansion_mul(_det6p, _val6p, _mul6p); + + /*---------------------------------- sum (-1)^P * VDi */ + INDX_TYPE + constexpr M1 = sub_alloc (N1, N2) ; + expansion _sum_1; + + INDX_TYPE + constexpr M2 = sub_alloc (N3, N4) ; + expansion _sum_2; + + INDX_TYPE + constexpr M3 = sub_alloc (N5, N6) ; + expansion _sum_3; + + if (_pivot % 2 == +0) + { + expansion_sub(_mul2p, _mul1p, _sum_1); + expansion_sub(_mul4p, _mul3p, _sum_2); + expansion_sub(_mul6p, _mul5p, _sum_3); + } + else + { + expansion_sub(_mul1p, _mul2p, _sum_1); + expansion_sub(_mul3p, _mul4p, _sum_2); + expansion_sub(_mul5p, _mul6p, _sum_3); + } + + expansion_add(_sum_1, _sum_2, _sum_3, _final); + } + + /*--------------------- "unitary" case, with Vi = +1. */ + + template < + size_t NA, size_t NB, size_t NC, + size_t ND, size_t NE, size_t NF, + size_t NG + > + __inline_call void unitary_det_6x6 ( + expansion const& _det1p , + expansion const& _det2p , + expansion const& _det3p , + expansion const& _det4p , + expansion const& _det5p , + expansion const& _det6p , + expansion & _final , + INDX_TYPE _pivot + ) + { + INDX_TYPE + constexpr N1 = sub_alloc (NA, NB) ; + expansion _sum_1; + + INDX_TYPE + constexpr N2 = sub_alloc (NC, ND) ; + expansion _sum_2; + + INDX_TYPE + constexpr M3 = sub_alloc (NE, NF) ; + expansion _sum_3; + + if (_pivot % 2 == +0) + { + expansion_sub(_det2p, _det1p, _sum_1); + expansion_sub(_det4p, _det3p, _sum_2); + expansion_sub(_det6p, _det5p, _sum_3); + } + else + { + expansion_sub(_det1p, _det2p, _sum_1); + expansion_sub(_det3p, _det4p, _sum_2); + expansion_sub(_det5p, _det6p, _sum_3); + } + + expansion_add(_sum_1, _sum_2, _sum_3, _final); + } + # undef REAL_TYPE # undef INDX_TYPE diff --git a/geompred.hpp b/geompred.hpp index 368d747..e4b6bea 100644 --- a/geompred.hpp +++ b/geompred.hpp @@ -31,11 +31,11 @@ * ------------------------------------------------------------ * - * Last updated: 01 March, 2020 + * Last updated: 11 May., 2024 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * ------------------------------------------------------------ diff --git a/mpfloats.hpp b/mpfloats.hpp index 21ffa94..46d8a21 100644 --- a/mpfloats.hpp +++ b/mpfloats.hpp @@ -31,11 +31,11 @@ * ------------------------------------------------------------ * - * Last updated: 11 April, 2020 + * Last updated: 11 May., 2024 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * ------------------------------------------------------------ diff --git a/predicate/bisect_k.hpp b/predicate/bisect_k.hpp index a3c32bf..bb3a730 100644 --- a/predicate/bisect_k.hpp +++ b/predicate/bisect_k.hpp @@ -22,20 +22,24 @@ * how they can obtain it for free, then you are not * required to make any arrangement with me.) * - * Disclaimer: Neither I nor: Columbia University, The - * Massachusetts Institute of Technology, The - * University of Sydney, nor The National Aeronautics - * and Space Administration warrant this code in any - * way whatsoever. This code is provided "as-is" to be - * used at your own risk. + * Disclaimer: Neither I nor THE CONTRIBUTORS warrant + * this code in any way whatsoever. This code is + * provided "as-is" to be used at your own risk. + * + * THE CONTRIBUTORS include: + * (a) The University of Sydney + * (b) The Massachusetts Institute of Technology + * (c) Columbia University + * (d) The National Aeronautics & Space Administration + * (e) Los Alamos National Laboratory * -------------------------------------------------------- * - * Last updated: 14 April, 2020 + * Last updated: 30 Apr., 2020 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * -------------------------------------------------------- @@ -624,5 +628,343 @@ return _sgn ; } + /* + -------------------------------------------------------- + * + * Compute an exact orientation wrt. bisector using + * multi-precision expansions, a'la shewchuk + * + * |c-a|**2 - wa = |c-b|**2 - wb + * + * This is the unweighted "bisect" predicate in E^4. + * + -------------------------------------------------------- + */ + + __normal_call REAL_TYPE bisect4d_e ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + bool_type &_OK + ) + { + /*--------------- bisect4d predicate, "exact" version */ + mp::expansion< 2 > _ab_xx_, _ab_yy_, + _ab_zz_, _ab_qq_, + _ac_xx_, _ac_yy_, + _ac_zz_, _ac_qq_, + _bc_xx_, _bc_yy_, + _bc_zz_, _bc_qq_; + mp::expansion< 4 > _tt_xx_, _tt_yy_, + _tt_zz_, _tt_qq_; + mp::expansion< 64> _absum_; + + _OK = true; + + /*----------------------------------- compute: d(p,q) */ + _ab_xx_.from_sub(_pa[0], _pb[0]); + _ab_yy_.from_sub(_pa[1], _pb[1]); + _ab_zz_.from_sub(_pa[2], _pb[2]); + _ab_qq_.from_sub(_pa[3], _pb[3]); + + _ac_xx_.from_sub(_pa[0], _pc[0]); + _ac_yy_.from_sub(_pa[1], _pc[1]); + _ac_zz_.from_sub(_pa[2], _pc[2]); + _ac_qq_.from_sub(_pa[3], _pc[3]); + + _bc_xx_.from_sub(_pb[0], _pc[0]); + _bc_yy_.from_sub(_pb[1], _pc[1]); + _bc_zz_.from_sub(_pb[2], _pc[2]); + _bc_qq_.from_sub(_pb[3], _pc[3]); + + mp::expansion_add(_ac_xx_, _bc_xx_, + _tt_xx_) ; + mp::expansion_add(_ac_yy_, _bc_yy_, + _tt_yy_) ; + mp::expansion_add(_ac_zz_, _bc_zz_, + _tt_zz_) ; + mp::expansion_add(_ac_qq_, _bc_qq_, + _tt_qq_) ; + + mp::expansion_dot(_ab_xx_, _tt_xx_, + _ab_yy_, _tt_yy_, + _ab_zz_, _tt_zz_, + _ab_qq_, _tt_qq_, + _absum_) ; + + /*----------------------------------- return signed d */ + return mp::expansion_est(_absum_) ; + } + + __normal_call REAL_TYPE bisect4d_i ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + bool_type &_OK + ) + { + /*--------------- bisect4d predicate, "bound" version */ + ia_flt _acx, _acy, _acz, _acq , + _bcx, _bcy, _bcz, _bcq , + _abx, _aby, _abz, _abq ; + ia_flt _sgn; + + ia_rnd _rnd; // up rounding! + + _abx.from_sub(_pa[0], _pb[0]) ; // coord. diff. + _aby.from_sub(_pa[1], _pb[1]) ; + _abz.from_sub(_pa[2], _pb[2]) ; + _abq.from_sub(_pa[3], _pb[3]) ; + + _acx.from_sub(_pa[0], _pc[0]) ; + _acy.from_sub(_pa[1], _pc[1]) ; + _acz.from_sub(_pa[2], _pc[2]) ; + _acq.from_sub(_pa[3], _pc[3]) ; + + _bcx.from_sub(_pb[0], _pc[0]) ; + _bcy.from_sub(_pb[1], _pc[1]) ; + _bcz.from_sub(_pb[2], _pc[2]) ; + _bcq.from_sub(_pb[3], _pc[3]) ; + + _sgn = (_abx * (_acx + _bcx)) + + (_aby * (_acy + _bcy)) + + (_abz * (_acz + _bcz)) + + (_abq * (_acq + _bcq)) ; + + _OK = + _sgn.lo() >= (REAL_TYPE)0. + || _sgn.up() <= (REAL_TYPE)0. ; + + return ( _sgn.mid () ) ; + } + + __normal_call REAL_TYPE bisect4d_f ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + bool_type &_OK + ) + { + /*--------------- bisect4d predicate, "float" version */ + REAL_TYPE static const _ER = + + 7. * std::pow(mp::_epsilon, 1) ; + + REAL_TYPE _acx, _acy, _acz, _acq; + REAL_TYPE _bcx, _bcy, _bcz, _bcq; + REAL_TYPE _acsqr, _bcsqr ; + + REAL_TYPE _sgn, _FT ; + + REAL_TYPE _ACSQR, _BCSQR ; + + _acx = _pa [0] - _pc [0] ; // coord. diff. + _acy = _pa [1] - _pc [1] ; + _acz = _pa [2] - _pc [2] ; + _acq = _pa [3] - _pc [3] ; + + _bcx = _pb [0] - _pc [0] ; + _bcy = _pb [1] - _pc [1] ; + _bcz = _pb [2] - _pc [2] ; + _bcq = _pb [3] - _pc [3] ; + + _acsqr = _acx * _acx + + _acy * _acy + + _acz * _acz + + _acq * _acq ; + _bcsqr = _bcx * _bcx + + _bcy * _bcy + + _bcz * _bcz + + _bcq * _bcq ; + + _ACSQR = std::abs(_acsqr); + _BCSQR = std::abs(_bcsqr); + + _FT = _ACSQR + _BCSQR ; // roundoff tol + _FT *= _ER ; + + _sgn = _acsqr - _bcsqr ; // d_ab - d_bc + + _OK = _sgn > +_FT || _sgn < -_FT ; + + return _sgn ; + } + + /* + -------------------------------------------------------- + * + * Compute an exact orientation wrt. bisector using + * multi-precision expansions, a'la shewchuk + * + * |c-a|**2 - wa = |c-b|**2 - wb + * + * This is the weighted "bisect" predicate in E^4. + * + -------------------------------------------------------- + */ + + __normal_call REAL_TYPE bisect4w_e ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + bool_type &_OK + ) + { + /*--------------- bisect4w predicate, "exact" version */ + mp::expansion< 2 > _ab_xx_, _ab_yy_, + _ab_zz_, _ab_qq_, + _ab_ww_, + _ac_xx_, _ac_yy_, + _ac_zz_, _ac_qq_, + _bc_xx_, _bc_yy_, + _bc_zz_, _bc_qq_; + mp::expansion< 4 > _tt_xx_, _tt_yy_, + _tt_zz_, _tt_qq_; + mp::expansion< 64> _ttsum_; + mp::expansion< 66> _absum_; + + _OK = true; + + /*----------------------------------- compute: d(p,q) */ + _ab_xx_.from_sub(_pa[0], _pb[0]); + _ab_yy_.from_sub(_pa[1], _pb[1]); + _ab_zz_.from_sub(_pa[2], _pb[2]); + _ab_qq_.from_sub(_pa[3], _pb[3]); + _ab_ww_.from_sub(_pa[4], _pb[4]); + + _ac_xx_.from_sub(_pa[0], _pc[0]); + _ac_yy_.from_sub(_pa[1], _pc[1]); + _ac_zz_.from_sub(_pa[2], _pc[2]); + _ac_qq_.from_sub(_pa[3], _pc[3]); + + _bc_xx_.from_sub(_pb[0], _pc[0]); + _bc_yy_.from_sub(_pb[1], _pc[1]); + _bc_zz_.from_sub(_pb[2], _pc[2]); + _bc_qq_.from_sub(_pb[3], _pc[3]); + + mp::expansion_add(_ac_xx_, _bc_xx_, + _tt_xx_) ; + mp::expansion_add(_ac_yy_, _bc_yy_, + _tt_yy_) ; + mp::expansion_add(_ac_zz_, _bc_zz_, + _tt_zz_) ; + mp::expansion_add(_ac_qq_, _bc_qq_, + _tt_qq_) ; + + mp::expansion_dot(_ab_xx_, _tt_xx_, + _ab_yy_, _tt_yy_, + _ab_zz_, _tt_zz_, + _ab_qq_, _tt_qq_, + _ttsum_) ; + + mp::expansion_sub(_ttsum_, _ab_ww_, + _absum_) ; + + /*----------------------------------- return signed d */ + return mp::expansion_est(_absum_) ; + } + + __normal_call REAL_TYPE bisect4w_i ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + bool_type &_OK + ) + { + /*--------------- bisect4w predicate, "bound" version */ + ia_flt _acx, _acy, _acz, _acq, + _abw, + _bcx, _bcy, _bcz, _bcq, + _abx, _aby, _abz, _abq; + ia_flt _sgn; + + ia_rnd _rnd; // up rounding! + + _abx.from_sub(_pa[0], _pb[0]) ; // coord. diff. + _aby.from_sub(_pa[1], _pb[1]) ; + _abz.from_sub(_pa[2], _pb[2]) ; + _abq.from_sub(_pa[3], _pb[3]) ; + _abw.from_sub(_pa[4], _pb[4]) ; + + _acx.from_sub(_pa[0], _pc[0]) ; + _acy.from_sub(_pa[1], _pc[1]) ; + _acz.from_sub(_pa[2], _pc[2]) ; + _acq.from_sub(_pa[3], _pc[3]) ; + + _bcx.from_sub(_pb[0], _pc[0]) ; + _bcy.from_sub(_pb[1], _pc[1]) ; + _bcz.from_sub(_pb[2], _pc[2]) ; + _bcq.from_sub(_pb[3], _pc[3]) ; + + _sgn = (_abx * (_acx + _bcx)) + + (_aby * (_acy + _bcy)) + + (_abz * (_acz + _bcz)) + + (_abq * (_acq + _bcq)) ; + + _sgn-= _abw ; + + _OK = + _sgn.lo() >= (REAL_TYPE)0. + || _sgn.up() <= (REAL_TYPE)0. ; + + return ( _sgn.mid () ) ; + } + + __normal_call REAL_TYPE bisect4w_f ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + bool_type &_OK + ) + { + /*--------------- bisect4w predicate, "float" version */ + REAL_TYPE static const _ER = + + 8. * std::pow(mp::_epsilon, 1) ; + + REAL_TYPE _acx, _acy, _acz, _acq; + REAL_TYPE _bcx, _bcy, _bcz, _bcq; + REAL_TYPE _acsqr, _bcsqr ; + REAL_TYPE _a_sum, _b_sum ; + + REAL_TYPE _A_SUM, _B_SUM ; + + REAL_TYPE _sgn, _FT; + + _acx = _pa [0] - _pc [0] ; // coord. diff. + _acy = _pa [1] - _pc [1] ; + _acz = _pa [2] - _pc [2] ; + _acq = _pa [3] - _pc [3] ; + + _bcx = _pb [0] - _pc [0] ; + _bcy = _pb [1] - _pc [1] ; + _bcz = _pb [2] - _pc [2] ; + _bcq = _pb [3] - _pc [3] ; + + _acsqr = _acx * _acx + + _acy * _acy + + _acz * _acz + + _acq * _acq ; + _bcsqr = _bcx * _bcx + + _bcy * _bcy + + _bcz * _bcz + + _bcq * _bcq ; + + _a_sum = _acsqr - _pa[4] ; + _b_sum = _bcsqr - _pb[4] ; + + _A_SUM = std::abs(_acsqr) + + std::abs(_pa[4]); + _B_SUM = std::abs(_bcsqr) + + std::abs(_pb[4]); + + _FT = _A_SUM + _B_SUM ; // roundoff tol + _FT *= _ER ; + + _sgn = _a_sum - _b_sum ; // d_ab - d_bc + + _OK = _sgn > +_FT || _sgn < -_FT ; + + return _sgn ; + } + diff --git a/predicate/inball_k.hpp b/predicate/inball_k.hpp index 5f716fd..12a85e3 100644 --- a/predicate/inball_k.hpp +++ b/predicate/inball_k.hpp @@ -22,20 +22,24 @@ * how they can obtain it for free, then you are not * required to make any arrangement with me.) * - * Disclaimer: Neither I nor: Columbia University, The - * Massachusetts Institute of Technology, The - * University of Sydney, nor The National Aeronautics - * and Space Administration warrant this code in any - * way whatsoever. This code is provided "as-is" to be - * used at your own risk. + * Disclaimer: Neither I nor THE CONTRIBUTORS warrant + * this code in any way whatsoever. This code is + * provided "as-is" to be used at your own risk. + * + * THE CONTRIBUTORS include: + * (a) The University of Sydney + * (b) The Massachusetts Institute of Technology + * (c) Columbia University + * (d) The National Aeronautics & Space Administration + * (e) Los Alamos National Laboratory * -------------------------------------------------------- * - * Last updated: 14 April, 2020 + * Last updated: 11 May., 2024 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * -------------------------------------------------------- @@ -1566,5 +1570,1009 @@ return ( _d44 ) ; } + /* + -------------------------------------------------------- + * + * Compute an exact determinant using multi-precision + * expansions, a'la shewchuk + * + * | ax ay az aq dot(a, a) +1. | + * | bx by bz bq dot(b, b) +1. | + * | cx cy cz cq dot(c, c) +1. | + * | dx dy dz dq dot(d, d) +1. | + * | ex ey ez eq dot(e, e) +1. | + * | fx fy fz fq dot(f, f) +1. | + * + * This is the unweighted "in-ball" predicate in E^4. + * + -------------------------------------------------------- + */ + + __normal_call REAL_TYPE inball4d_e ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + __const_ptr(REAL_TYPE) _pd , + __const_ptr(REAL_TYPE) _pe , + __const_ptr(REAL_TYPE) _pf , + bool_type &_OK + ) + { + /*--------------- inball4d predicate, "exact" version */ + mp::expansion< 8 > _a_lift, _b_lift, + _c_lift, _d_lift, + _e_lift, _f_lift; + mp::expansion< 4 > _d2_ab_, _d2_ac_, + _d2_ad_, _d2_ae_, + _d2_af_, + _d2_bc_, _d2_bd_, + _d2_be_, _d2_bf_, + _d2_cd_, _d2_ce_, + _d2_cf_, + _d2_de_, _d2_df_, + _d2_ef_; + mp::expansion< 24> _d3_abc, _d3_abd, + _d3_abe, _d3_abf, + _d3_acd, _d3_ace, + _d3_acf, + _d3_ade, _d3_adf, + _d3_aef, + _d3_bcd, _d3_bce, + _d3_bcf, + _d3_bde, _d3_bdf, + _d3_bef, + _d3_cde, _d3_cdf, + _d3_cef, _d3_def; + mp::expansion<192> _d4abcd, _d4abce, + _d4abcf, + _d4abde, _d4abdf, + _d4abef, + _d4acde, _d4acdf, + _d4acef, _d4adef, + _d4bcde, _d4bcdf, + _d4bcef, _d4bdef, + _d4cdef; + mp::expansion<960> _5bcdef, _5acdef, + _5abdef, _5abcef, + _5abcdf, _5abcde; + // try not to blow the stack... + auto _d6full = new mp::expansion<92160>() ; + + _OK = true; + + mp::expansion< 1 > _pa_zz_(_pa[ 2]); + mp::expansion< 1 > _pb_zz_(_pb[ 2]); + mp::expansion< 1 > _pc_zz_(_pc[ 2]); + mp::expansion< 1 > _pd_zz_(_pd[ 2]); + mp::expansion< 1 > _pe_zz_(_pe[ 2]); + mp::expansion< 1 > _pf_zz_(_pf[ 2]); + + mp::expansion< 1 > _pa_qq_(_pa[ 3]); + mp::expansion< 1 > _pb_qq_(_pb[ 3]); + mp::expansion< 1 > _pc_qq_(_pc[ 3]); + mp::expansion< 1 > _pd_qq_(_pd[ 3]); + mp::expansion< 1 > _pe_qq_(_pe[ 3]); + mp::expansion< 1 > _pf_qq_(_pf[ 3]); + + /*-------------------------------------- lifted terms */ + mp::expansion_add( + mp::expansion_from_sqr(_pa[ 0]), + mp::expansion_from_sqr(_pa[ 1]), + mp::expansion_from_sqr(_pa[ 2]), + mp::expansion_from_sqr(_pa[ 3]), + _a_lift ) ; + + mp::expansion_add( + mp::expansion_from_sqr(_pb[ 0]), + mp::expansion_from_sqr(_pb[ 1]), + mp::expansion_from_sqr(_pb[ 2]), + mp::expansion_from_sqr(_pb[ 3]), + _b_lift ) ; + + mp::expansion_add( + mp::expansion_from_sqr(_pc[ 0]), + mp::expansion_from_sqr(_pc[ 1]), + mp::expansion_from_sqr(_pc[ 2]), + mp::expansion_from_sqr(_pc[ 3]), + _c_lift ) ; + + mp::expansion_add( + mp::expansion_from_sqr(_pd[ 0]), + mp::expansion_from_sqr(_pd[ 1]), + mp::expansion_from_sqr(_pd[ 2]), + mp::expansion_from_sqr(_pd[ 3]), + _d_lift ) ; + + mp::expansion_add( + mp::expansion_from_sqr(_pe[ 0]), + mp::expansion_from_sqr(_pe[ 1]), + mp::expansion_from_sqr(_pe[ 2]), + mp::expansion_from_sqr(_pe[ 3]), + _e_lift ) ; + + mp::expansion_add( + mp::expansion_from_sqr(_pf[ 0]), + mp::expansion_from_sqr(_pf[ 1]), + mp::expansion_from_sqr(_pf[ 2]), + mp::expansion_from_sqr(_pf[ 3]), + _f_lift ) ; + + /*-------------------------------------- 2 x 2 minors */ + compute_det_2x2(_pa[ 0], _pa[ 1], + _pb[ 0], _pb[ 1], + _d2_ab_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pc[ 0], _pc[ 1], + _d2_ac_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pd[ 0], _pd[ 1], + _d2_ad_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pe[ 0], _pe[ 1], + _d2_ae_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pf[ 0], _pf[ 1], + _d2_af_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pc[ 0], _pc[ 1], + _d2_bc_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pd[ 0], _pd[ 1], + _d2_bd_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pe[ 0], _pe[ 1], + _d2_be_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pf[ 0], _pf[ 1], + _d2_bf_ ) ; + + compute_det_2x2(_pc[ 0], _pc[ 1], + _pd[ 0], _pd[ 1], + _d2_cd_ ) ; + + compute_det_2x2(_pc[ 0], _pc[ 1], + _pe[ 0], _pe[ 1], + _d2_ce_ ) ; + + compute_det_2x2(_pc[ 0], _pc[ 1], + _pf[ 0], _pf[ 1], + _d2_cf_ ) ; + + compute_det_2x2(_pd[ 0], _pd[ 1], + _pe[ 0], _pe[ 1], + _d2_de_ ) ; + + compute_det_2x2(_pd[ 0], _pd[ 1], + _pf[ 0], _pf[ 1], + _d2_df_ ) ; + + compute_det_2x2(_pe[ 0], _pe[ 1], + _pf[ 0], _pf[ 1], + _d2_ef_ ) ; + + /*-------------------------------------- 3 x 3 minors */ + compute_det_3x3(_d2_bc_, _pa_zz_, + _d2_ac_, _pb_zz_, + _d2_ab_, _pc_zz_, + _d3_abc, +3) ; + + compute_det_3x3(_d2_bd_, _pa_zz_, + _d2_ad_, _pb_zz_, + _d2_ab_, _pd_zz_, + _d3_abd, +3) ; + + compute_det_3x3(_d2_be_, _pa_zz_, + _d2_ae_, _pb_zz_, + _d2_ab_, _pe_zz_, + _d3_abe, +3) ; + + compute_det_3x3(_d2_bf_, _pa_zz_, + _d2_af_, _pb_zz_, + _d2_ab_, _pf_zz_, + _d3_abf, +3) ; + + compute_det_3x3(_d2_cd_, _pa_zz_, + _d2_ad_, _pc_zz_, + _d2_ac_, _pd_zz_, + _d3_acd, +3) ; + + compute_det_3x3(_d2_ce_, _pa_zz_, + _d2_ae_, _pc_zz_, + _d2_ac_, _pe_zz_, + _d3_ace, +3) ; + + compute_det_3x3(_d2_cf_, _pa_zz_, + _d2_af_, _pc_zz_, + _d2_ac_, _pf_zz_, + _d3_acf, +3) ; + + compute_det_3x3(_d2_de_, _pa_zz_, + _d2_ae_, _pd_zz_, + _d2_ad_, _pe_zz_, + _d3_ade, +3) ; + + compute_det_3x3(_d2_df_, _pa_zz_, + _d2_af_, _pd_zz_, + _d2_ad_, _pf_zz_, + _d3_adf, +3) ; + + compute_det_3x3(_d2_ef_, _pa_zz_, + _d2_af_, _pe_zz_, + _d2_ae_, _pf_zz_, + _d3_aef, +3) ; + + compute_det_3x3(_d2_cd_, _pb_zz_, + _d2_bd_, _pc_zz_, + _d2_bc_, _pd_zz_, + _d3_bcd, +3) ; + + compute_det_3x3(_d2_ce_, _pb_zz_, + _d2_be_, _pc_zz_, + _d2_bc_, _pe_zz_, + _d3_bce, +3) ; + + compute_det_3x3(_d2_cf_, _pb_zz_, + _d2_bf_, _pc_zz_, + _d2_bc_, _pf_zz_, + _d3_bcf, +3) ; + + compute_det_3x3(_d2_de_, _pb_zz_, + _d2_be_, _pd_zz_, + _d2_bd_, _pe_zz_, + _d3_bde, +3) ; + + compute_det_3x3(_d2_df_, _pb_zz_, + _d2_bf_, _pd_zz_, + _d2_bd_, _pf_zz_, + _d3_bdf, +3) ; + + compute_det_3x3(_d2_ef_, _pb_zz_, + _d2_bf_, _pe_zz_, + _d2_be_, _pf_zz_, + _d3_bef, +3) ; + + compute_det_3x3(_d2_de_, _pc_zz_, + _d2_ce_, _pd_zz_, + _d2_cd_, _pe_zz_, + _d3_cde, +3) ; + + compute_det_3x3(_d2_df_, _pc_zz_, + _d2_cf_, _pd_zz_, + _d2_cd_, _pf_zz_, + _d3_cdf, +3) ; + + compute_det_3x3(_d2_ef_, _pc_zz_, + _d2_cf_, _pe_zz_, + _d2_ce_, _pf_zz_, + _d3_cef, +3) ; + + compute_det_3x3(_d2_ef_, _pd_zz_, + _d2_df_, _pe_zz_, + _d2_de_, _pf_zz_, + _d3_def, +3) ; + + /*-------------------------------------- 4 x 4 minors */ + compute_det_4x4(_d3_bcd, _pa_qq_, + _d3_acd, _pb_qq_, + _d3_abd, _pc_qq_, + _d3_abc, _pd_qq_, + _d4abcd, +4) ; + + compute_det_4x4(_d3_bce, _pa_qq_, + _d3_ace, _pb_qq_, + _d3_abe, _pc_qq_, + _d3_abc, _pe_qq_, + _d4abce, +4) ; + + compute_det_4x4(_d3_bcf, _pa_qq_, + _d3_acf, _pb_qq_, + _d3_abf, _pc_qq_, + _d3_abc, _pf_qq_, + _d4abcf, +4) ; + + compute_det_4x4(_d3_bde, _pa_qq_, + _d3_ade, _pb_qq_, + _d3_abe, _pd_qq_, + _d3_abd, _pe_qq_, + _d4abde, +4) ; + + compute_det_4x4(_d3_bdf, _pa_qq_, + _d3_adf, _pb_qq_, + _d3_abf, _pd_qq_, + _d3_abd, _pf_qq_, + _d4abdf, +4) ; + + compute_det_4x4(_d3_bef, _pa_qq_, + _d3_aef, _pb_qq_, + _d3_abf, _pe_qq_, + _d3_abe, _pf_qq_, + _d4abef, +4) ; + + compute_det_4x4(_d3_cde, _pa_qq_, + _d3_ade, _pc_qq_, + _d3_ace, _pd_qq_, + _d3_acd, _pe_qq_, + _d4acde, +4) ; + + compute_det_4x4(_d3_cdf, _pa_qq_, + _d3_adf, _pc_qq_, + _d3_acf, _pd_qq_, + _d3_acd, _pf_qq_, + _d4acdf, +4) ; + + compute_det_4x4(_d3_cef, _pa_qq_, + _d3_aef, _pc_qq_, + _d3_acf, _pe_qq_, + _d3_ace, _pf_qq_, + _d4acef, +4) ; + + compute_det_4x4(_d3_def, _pa_qq_, + _d3_aef, _pd_qq_, + _d3_adf, _pe_qq_, + _d3_ade, _pf_qq_, + _d4adef, +4) ; + + compute_det_4x4(_d3_cde, _pb_qq_, + _d3_bde, _pc_qq_, + _d3_bce, _pd_qq_, + _d3_bcd, _pe_qq_, + _d4bcde, +4) ; + + compute_det_4x4(_d3_cdf, _pb_qq_, + _d3_bdf, _pc_qq_, + _d3_bcf, _pd_qq_, + _d3_bcd, _pf_qq_, + _d4bcdf, +4) ; + + compute_det_4x4(_d3_cef, _pb_qq_, + _d3_bef, _pc_qq_, + _d3_bcf, _pe_qq_, + _d3_bce, _pf_qq_, + _d4bcef, +4) ; + + compute_det_4x4(_d3_def, _pb_qq_, + _d3_bef, _pd_qq_, + _d3_bdf, _pe_qq_, + _d3_bde, _pf_qq_, + _d4bdef, +4) ; + + compute_det_4x4(_d3_def, _pc_qq_, + _d3_cef, _pd_qq_, + _d3_cdf, _pe_qq_, + _d3_cde, _pf_qq_, + _d4cdef, +4) ; + + /*-------------------------------------- 5 x 5 minors */ + unitary_det_5x5(_d4bcde, _d4acde, + _d4abde, _d4abce, + _d4abcd, + _5abcde, +5) ; + + unitary_det_5x5(_d4bcdf, _d4acdf, + _d4abdf, _d4abcf, + _d4abcd, + _5abcdf, +5) ; + + unitary_det_5x5(_d4bcef, _d4acef, + _d4abef, _d4abcf, + _d4abce, + _5abcef, +5) ; + + unitary_det_5x5(_d4bdef, _d4adef, + _d4abef, _d4abdf, + _d4abde, + _5abdef, +5) ; + + unitary_det_5x5(_d4cdef, _d4adef, + _d4acef, _d4acdf, + _d4acde, + _5acdef, +5) ; + + unitary_det_5x5(_d4cdef, _d4bdef, + _d4bcef, _d4bcdf, + _d4bcde, + _5bcdef, +5) ; + + /*-------------------------------------- 6 x 6 result */ + compute_det_6x6(_5bcdef, _a_lift, + _5acdef, _b_lift, + _5abdef, _c_lift, + _5abcef, _d_lift, + _5abcdf, _e_lift, + _5abcde, _f_lift, + *_d6full, +5) ; + + /*-------------------------------------- leading det. */ + REAL_TYPE _d66 = mp::expansion_est(*_d6full) ; + + delete _d6full ; return _d66 ; + } + + __normal_call REAL_TYPE inball4d_i ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + __const_ptr(REAL_TYPE) _pd , + __const_ptr(REAL_TYPE) _pe , + __const_ptr(REAL_TYPE) _pf , + bool_type &_OK + ) + { + /*--------------- inball4d predicate, "bound" version */ + ia_flt _afx, _afy, _afz , + _afq, _ali, + _bfx, _bfy, _bfz , + _bfq, _bli, + _cfx, _cfy, _cfz , + _cfq, _cli, + _dfx, _dfy, _dfz , + _dfq, _dli, + _efx, _efy, _efz , + _efq, _eli; + + ia_flt _afxbfy, _bfxafy , + _afxcfy, _cfxafy , + _bfxcfy, _cfxbfy , + _cfxdfy, _dfxcfy , + _dfxafy, _afxdfy , + _bfxdfy, _dfxbfy ; + + ia_flt _ab_, _bc_, _cd_, _da_, + _ac_, _bd_; + + + ia_flt _abc, _bcd, _cda, _dab; + + + ia_flt _d55; + + ia_rnd _rnd; // up rounding! + + _afx.from_sub(_pa[0], _pf[0]) ; // coord. diff. + _afy.from_sub(_pa[1], _pf[1]) ; + _afz.from_sub(_pa[2], _pf[2]) ; + _afq.from_sub(_pa[3], _pf[3]) ; + + _bfx.from_sub(_pb[0], _pf[0]) ; + _bfy.from_sub(_pb[1], _pf[1]) ; + _bfz.from_sub(_pb[2], _pf[2]) ; + _bfq.from_sub(_pb[3], _pf[3]) ; + + _cfx.from_sub(_pc[0], _pf[0]) ; + _cfy.from_sub(_pc[1], _pf[1]) ; + _cfz.from_sub(_pc[2], _pf[2]) ; + _cfq.from_sub(_pc[3], _pf[3]) ; + + _dfx.from_sub(_pd[0], _pf[0]) ; + _dfy.from_sub(_pd[1], _pf[1]) ; + _dfz.from_sub(_pd[2], _pf[2]) ; + _dfq.from_sub(_pd[3], _pf[3]) ; + + _efx.from_sub(_pe[0], _pf[0]) ; + _efy.from_sub(_pe[1], _pf[1]) ; + _efz.from_sub(_pe[2], _pf[2]) ; + _efq.from_sub(_pe[3], _pf[3]) ; + + _ali = sqr (_afx) + sqr (_afy) // lifted terms + + sqr (_afz) + sqr (_afq) ; + + _bli = sqr (_bfx) + sqr (_bfy) + + sqr (_bfz) + sqr (_bfq) ; + + _cli = sqr (_cfx) + sqr (_cfy) + + sqr (_cfz) + sqr (_cfq) ; + + _dli = sqr (_dfx) + sqr (_dfy) + + sqr (_dfz) + sqr (_dfq) ; + + _eli = sqr (_efx) + sqr (_efy) + + sqr (_efz) + sqr (_efq) ; + + + /* + _aexbey = _aex * _bey ; // 2 x 2 minors + _bexaey = _bex * _aey ; + _ab_ = _aexbey - _bexaey ; + + _bexcey = _bex * _cey; + _cexbey = _cex * _bey; + _bc_ = _bexcey - _cexbey ; + + _cexdey = _cex * _dey; + _dexcey = _dex * _cey; + _cd_ = _cexdey - _dexcey ; + + _dexaey = _dex * _aey; + _aexdey = _aex * _dey; + _da_ = _dexaey - _aexdey ; + + _aexcey = _aex * _cey; + _cexaey = _cex * _aey; + _ac_ = _aexcey - _cexaey ; + + _bexdey = _bex * _dey; + _dexbey = _dex * _bey; + _bd_ = _bexdey - _dexbey ; + + + _abc = // 3 x 3 minors + _aez * _bc_ - _bez * _ac_ + + _cez * _ab_ ; + + _bcd = + _bez * _cd_ - _cez * _bd_ + + _dez * _bc_ ; + + _cda = + _cez * _da_ + _dez * _ac_ + + _aez * _cd_ ; + + _dab = + _dez * _ab_ + _aez * _bd_ + + _bez * _da_ ; + + + _d44 = // 5 x 5 result + _eli * _abcd + + _dli * _abce + - _cli * _deab + + _bli * _cdea + - _ali * _bcde ; + */ + + _OK = + _d55.lo() >= (REAL_TYPE)0. + ||_d55.up() <= (REAL_TYPE)0.; + + return ( _d55.mid() ) ; + } + + /* + -------------------------------------------------------- + * + * Compute an exact determinant using multi-precision + * expansions, a'la shewchuk + * + * | ax ay az aq dot(a, a) - aw +1. | + * | bx by bz bq dot(b, b) - bw +1. | + * | cx cy cz cq dot(c, c) - cw +1. | + * | dx dy dz dq dot(d, d) - dw +1. | + * | ex ey ez eq dot(e, e) - ew +1. | + * | fx fy fz fq dot(f, f) - fw +1. | + * + * This is the weighted "in-ball" predicate in E^4. + * + -------------------------------------------------------- + */ + + __normal_call REAL_TYPE inball4w_e ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + __const_ptr(REAL_TYPE) _pd , + __const_ptr(REAL_TYPE) _pe , + __const_ptr(REAL_TYPE) _pf , + bool_type &_OK + ) + { + /*--------------- inball4w predicate, "exact" version */ + mp::expansion< 9 > _a_lift, _b_lift, + _c_lift, _d_lift, + _e_lift, _f_lift; + mp::expansion< 8 > _t_lift; + mp::expansion< 4 > _d2_ab_, _d2_ac_, + _d2_ad_, _d2_ae_, + _d2_af_, + _d2_bc_, _d2_bd_, + _d2_be_, _d2_bf_, + _d2_cd_, _d2_ce_, + _d2_cf_, + _d2_de_, _d2_df_, + _d2_ef_; + mp::expansion< 24> _d3_abc, _d3_abd, + _d3_abe, _d3_abf, + _d3_acd, _d3_ace, + _d3_acf, + _d3_ade, _d3_adf, + _d3_aef, + _d3_bcd, _d3_bce, + _d3_bcf, + _d3_bde, _d3_bdf, + _d3_bef, + _d3_cde, _d3_cdf, + _d3_cef, _d3_def; + mp::expansion<192> _d4abcd, _d4abce, + _d4abcf, + _d4abde, _d4abdf, + _d4abef, + _d4acde, _d4acdf, + _d4acef, _d4adef, + _d4bcde, _d4bcdf, + _d4bcef, _d4bdef, + _d4cdef; + mp::expansion<960> _5bcdef, _5acdef, + _5abdef, _5abcef, + _5abcdf, _5abcde; + // try not to blow the stack... + auto _d6full = new mp::expansion<103680>(); + + _OK = true; + + mp::expansion< 1 > _pa_zz_(_pa[ 2]); + mp::expansion< 1 > _pb_zz_(_pb[ 2]); + mp::expansion< 1 > _pc_zz_(_pc[ 2]); + mp::expansion< 1 > _pd_zz_(_pd[ 2]); + mp::expansion< 1 > _pe_zz_(_pe[ 2]); + mp::expansion< 1 > _pf_zz_(_pf[ 2]); + + mp::expansion< 1 > _pa_qq_(_pa[ 3]); + mp::expansion< 1 > _pb_qq_(_pb[ 3]); + mp::expansion< 1 > _pc_qq_(_pc[ 3]); + mp::expansion< 1 > _pd_qq_(_pd[ 3]); + mp::expansion< 1 > _pe_qq_(_pe[ 3]); + mp::expansion< 1 > _pf_qq_(_pf[ 3]); + + /*-------------------------------------- lifted terms */ + mp::expansion_add( + mp::expansion_from_sqr(_pa[ 0]), + mp::expansion_from_sqr(_pa[ 1]), + mp::expansion_from_sqr(_pa[ 2]), + mp::expansion_from_sqr(_pa[ 3]), + _t_lift ) ; + mp::expansion_sub( + _t_lift , _pa[ 4] , _a_lift); + + mp::expansion_add( + mp::expansion_from_sqr(_pb[ 0]), + mp::expansion_from_sqr(_pb[ 1]), + mp::expansion_from_sqr(_pb[ 2]), + mp::expansion_from_sqr(_pb[ 3]), + _t_lift ) ; + mp::expansion_sub( + _t_lift , _pb[ 4] , _b_lift); + + mp::expansion_add( + mp::expansion_from_sqr(_pc[ 0]), + mp::expansion_from_sqr(_pc[ 1]), + mp::expansion_from_sqr(_pc[ 2]), + mp::expansion_from_sqr(_pc[ 3]), + _t_lift ) ; + mp::expansion_sub( + _t_lift , _pc[ 4] , _c_lift); + + mp::expansion_add( + mp::expansion_from_sqr(_pd[ 0]), + mp::expansion_from_sqr(_pd[ 1]), + mp::expansion_from_sqr(_pd[ 2]), + mp::expansion_from_sqr(_pd[ 3]), + _t_lift ) ; + mp::expansion_sub( + _t_lift , _pd[ 4] , _d_lift); + + mp::expansion_add( + mp::expansion_from_sqr(_pe[ 0]), + mp::expansion_from_sqr(_pe[ 1]), + mp::expansion_from_sqr(_pe[ 2]), + mp::expansion_from_sqr(_pe[ 3]), + _t_lift ) ; + mp::expansion_sub( + _t_lift , _pe[ 4] , _e_lift); + + mp::expansion_add( + mp::expansion_from_sqr(_pf[ 0]), + mp::expansion_from_sqr(_pf[ 1]), + mp::expansion_from_sqr(_pf[ 2]), + mp::expansion_from_sqr(_pf[ 3]), + _t_lift ) ; + mp::expansion_sub( + _t_lift , _pf[ 4] , _f_lift); + + /*-------------------------------------- 2 x 2 minors */ + compute_det_2x2(_pa[ 0], _pa[ 1], + _pb[ 0], _pb[ 1], + _d2_ab_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pc[ 0], _pc[ 1], + _d2_ac_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pd[ 0], _pd[ 1], + _d2_ad_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pe[ 0], _pe[ 1], + _d2_ae_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pf[ 0], _pf[ 1], + _d2_af_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pc[ 0], _pc[ 1], + _d2_bc_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pd[ 0], _pd[ 1], + _d2_bd_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pe[ 0], _pe[ 1], + _d2_be_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pf[ 0], _pf[ 1], + _d2_bf_ ) ; + + compute_det_2x2(_pc[ 0], _pc[ 1], + _pd[ 0], _pd[ 1], + _d2_cd_ ) ; + + compute_det_2x2(_pc[ 0], _pc[ 1], + _pe[ 0], _pe[ 1], + _d2_ce_ ) ; + + compute_det_2x2(_pc[ 0], _pc[ 1], + _pf[ 0], _pf[ 1], + _d2_cf_ ) ; + + compute_det_2x2(_pd[ 0], _pd[ 1], + _pe[ 0], _pe[ 1], + _d2_de_ ) ; + + compute_det_2x2(_pd[ 0], _pd[ 1], + _pf[ 0], _pf[ 1], + _d2_df_ ) ; + + compute_det_2x2(_pe[ 0], _pe[ 1], + _pf[ 0], _pf[ 1], + _d2_ef_ ) ; + + /*-------------------------------------- 3 x 3 minors */ + compute_det_3x3(_d2_bc_, _pa_zz_, + _d2_ac_, _pb_zz_, + _d2_ab_, _pc_zz_, + _d3_abc, +3) ; + + compute_det_3x3(_d2_bd_, _pa_zz_, + _d2_ad_, _pb_zz_, + _d2_ab_, _pd_zz_, + _d3_abd, +3) ; + + compute_det_3x3(_d2_be_, _pa_zz_, + _d2_ae_, _pb_zz_, + _d2_ab_, _pe_zz_, + _d3_abe, +3) ; + + compute_det_3x3(_d2_bf_, _pa_zz_, + _d2_af_, _pb_zz_, + _d2_ab_, _pf_zz_, + _d3_abf, +3) ; + + compute_det_3x3(_d2_cd_, _pa_zz_, + _d2_ad_, _pc_zz_, + _d2_ac_, _pd_zz_, + _d3_acd, +3) ; + + compute_det_3x3(_d2_ce_, _pa_zz_, + _d2_ae_, _pc_zz_, + _d2_ac_, _pe_zz_, + _d3_ace, +3) ; + + compute_det_3x3(_d2_cf_, _pa_zz_, + _d2_af_, _pc_zz_, + _d2_ac_, _pf_zz_, + _d3_acf, +3) ; + + compute_det_3x3(_d2_de_, _pa_zz_, + _d2_ae_, _pd_zz_, + _d2_ad_, _pe_zz_, + _d3_ade, +3) ; + + compute_det_3x3(_d2_df_, _pa_zz_, + _d2_af_, _pd_zz_, + _d2_ad_, _pf_zz_, + _d3_adf, +3) ; + + compute_det_3x3(_d2_ef_, _pa_zz_, + _d2_af_, _pe_zz_, + _d2_ae_, _pf_zz_, + _d3_aef, +3) ; + + compute_det_3x3(_d2_cd_, _pb_zz_, + _d2_bd_, _pc_zz_, + _d2_bc_, _pd_zz_, + _d3_bcd, +3) ; + + compute_det_3x3(_d2_ce_, _pb_zz_, + _d2_be_, _pc_zz_, + _d2_bc_, _pe_zz_, + _d3_bce, +3) ; + + compute_det_3x3(_d2_cf_, _pb_zz_, + _d2_bf_, _pc_zz_, + _d2_bc_, _pf_zz_, + _d3_bcf, +3) ; + + compute_det_3x3(_d2_de_, _pb_zz_, + _d2_be_, _pd_zz_, + _d2_bd_, _pe_zz_, + _d3_bde, +3) ; + + compute_det_3x3(_d2_df_, _pb_zz_, + _d2_bf_, _pd_zz_, + _d2_bd_, _pf_zz_, + _d3_bdf, +3) ; + + compute_det_3x3(_d2_ef_, _pb_zz_, + _d2_bf_, _pe_zz_, + _d2_be_, _pf_zz_, + _d3_bef, +3) ; + + compute_det_3x3(_d2_de_, _pc_zz_, + _d2_ce_, _pd_zz_, + _d2_cd_, _pe_zz_, + _d3_cde, +3) ; + + compute_det_3x3(_d2_df_, _pc_zz_, + _d2_cf_, _pd_zz_, + _d2_cd_, _pf_zz_, + _d3_cdf, +3) ; + + compute_det_3x3(_d2_ef_, _pc_zz_, + _d2_cf_, _pe_zz_, + _d2_ce_, _pf_zz_, + _d3_cef, +3) ; + + compute_det_3x3(_d2_ef_, _pd_zz_, + _d2_df_, _pe_zz_, + _d2_de_, _pf_zz_, + _d3_def, +3) ; + + /*-------------------------------------- 4 x 4 minors */ + compute_det_4x4(_d3_bcd, _pa_qq_, + _d3_acd, _pb_qq_, + _d3_abd, _pc_qq_, + _d3_abc, _pd_qq_, + _d4abcd, +4) ; + + compute_det_4x4(_d3_bce, _pa_qq_, + _d3_ace, _pb_qq_, + _d3_abe, _pc_qq_, + _d3_abc, _pe_qq_, + _d4abce, +4) ; + + compute_det_4x4(_d3_bcf, _pa_qq_, + _d3_acf, _pb_qq_, + _d3_abf, _pc_qq_, + _d3_abc, _pf_qq_, + _d4abcf, +4) ; + + compute_det_4x4(_d3_bde, _pa_qq_, + _d3_ade, _pb_qq_, + _d3_abe, _pd_qq_, + _d3_abd, _pe_qq_, + _d4abde, +4) ; + + compute_det_4x4(_d3_bdf, _pa_qq_, + _d3_adf, _pb_qq_, + _d3_abf, _pd_qq_, + _d3_abd, _pf_qq_, + _d4abdf, +4) ; + + compute_det_4x4(_d3_bef, _pa_qq_, + _d3_aef, _pb_qq_, + _d3_abf, _pe_qq_, + _d3_abe, _pf_qq_, + _d4abef, +4) ; + + compute_det_4x4(_d3_cde, _pa_qq_, + _d3_ade, _pc_qq_, + _d3_ace, _pd_qq_, + _d3_acd, _pe_qq_, + _d4acde, +4) ; + + compute_det_4x4(_d3_cdf, _pa_qq_, + _d3_adf, _pc_qq_, + _d3_acf, _pd_qq_, + _d3_acd, _pf_qq_, + _d4acdf, +4) ; + + compute_det_4x4(_d3_cef, _pa_qq_, + _d3_aef, _pc_qq_, + _d3_acf, _pe_qq_, + _d3_ace, _pf_qq_, + _d4acef, +4) ; + + compute_det_4x4(_d3_def, _pa_qq_, + _d3_aef, _pd_qq_, + _d3_adf, _pe_qq_, + _d3_ade, _pf_qq_, + _d4adef, +4) ; + + compute_det_4x4(_d3_cde, _pb_qq_, + _d3_bde, _pc_qq_, + _d3_bce, _pd_qq_, + _d3_bcd, _pe_qq_, + _d4bcde, +4) ; + + compute_det_4x4(_d3_cdf, _pb_qq_, + _d3_bdf, _pc_qq_, + _d3_bcf, _pd_qq_, + _d3_bcd, _pf_qq_, + _d4bcdf, +4) ; + + compute_det_4x4(_d3_cef, _pb_qq_, + _d3_bef, _pc_qq_, + _d3_bcf, _pe_qq_, + _d3_bce, _pf_qq_, + _d4bcef, +4) ; + + compute_det_4x4(_d3_def, _pb_qq_, + _d3_bef, _pd_qq_, + _d3_bdf, _pe_qq_, + _d3_bde, _pf_qq_, + _d4bdef, +4) ; + + compute_det_4x4(_d3_def, _pc_qq_, + _d3_cef, _pd_qq_, + _d3_cdf, _pe_qq_, + _d3_cde, _pf_qq_, + _d4cdef, +4) ; + + /*-------------------------------------- 5 x 5 minors */ + unitary_det_5x5(_d4bcde, _d4acde, + _d4abde, _d4abce, + _d4abcd, + _5abcde, +5) ; + + unitary_det_5x5(_d4bcdf, _d4acdf, + _d4abdf, _d4abcf, + _d4abcd, + _5abcdf, +5) ; + + unitary_det_5x5(_d4bcef, _d4acef, + _d4abef, _d4abcf, + _d4abce, + _5abcef, +5) ; + + unitary_det_5x5(_d4bdef, _d4adef, + _d4abef, _d4abdf, + _d4abde, + _5abdef, +5) ; + + unitary_det_5x5(_d4cdef, _d4adef, + _d4acef, _d4acdf, + _d4acde, + _5acdef, +5) ; + + unitary_det_5x5(_d4cdef, _d4bdef, + _d4bcef, _d4bcdf, + _d4bcde, + _5bcdef, +5) ; + + /*-------------------------------------- 6 x 6 result */ + compute_det_6x6(_5bcdef, _a_lift, + _5acdef, _b_lift, + _5abdef, _c_lift, + _5abcef, _d_lift, + _5abcdf, _e_lift, + _5abcde, _f_lift, + *_d6full, +5) ; + + /*-------------------------------------- leading det. */ + REAL_TYPE _d66 = mp::expansion_est(*_d6full) ; + + delete _d6full ; return _d66 ; + } + diff --git a/predicate/orient_k.hpp b/predicate/orient_k.hpp index 1c27883..fe6505d 100644 --- a/predicate/orient_k.hpp +++ b/predicate/orient_k.hpp @@ -22,20 +22,24 @@ * how they can obtain it for free, then you are not * required to make any arrangement with me.) * - * Disclaimer: Neither I nor: Columbia University, The - * Massachusetts Institute of Technology, The - * University of Sydney, nor The National Aeronautics - * and Space Administration warrant this code in any - * way whatsoever. This code is provided "as-is" to be - * used at your own risk. + * Disclaimer: Neither I nor THE CONTRIBUTORS warrant + * this code in any way whatsoever. This code is + * provided "as-is" to be used at your own risk. + * + * THE CONTRIBUTORS include: + * (a) The University of Sydney + * (b) The Massachusetts Institute of Technology + * (c) Columbia University + * (d) The National Aeronautics & Space Administration + * (e) Los Alamos National Laboratory * -------------------------------------------------------- * - * Last updated: 14 April, 2020 + * Last updated: 30 Apr., 2020 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * -------------------------------------------------------- @@ -391,5 +395,458 @@ return ( _sgn ) ; } + /* + -------------------------------------------------------- + * + * Compute an exact determinant using multi-precision + * expansions, a'la shewchuk + * + * | ax ay az aq +1. | + * | bx by bz bq +1. | + * | cx cy cz cq +1. | + * | dx dy dz dq +1. | + * | ex ey ez eq +1. | + * + * This is the planar "orientation" predicate in E^4. + * + -------------------------------------------------------- + */ + + __normal_call REAL_TYPE orient4d_e ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + __const_ptr(REAL_TYPE) _pd , + __const_ptr(REAL_TYPE) _pe , + bool_type &_OK + ) + { + /*--------------- orient4d predicate, "exact" version */ + mp::expansion< 4 > _d2_ab_, _d2_ac_, + _d2_ad_, _d2_ae_, + _d2_bc_, _d2_bd_, + _d2_be_, + _d2_cd_, _d2_ce_, + _d2_de_; + mp::expansion< 24> _d3_abc, _d3_abd, + _d3_abe, + _d3_acd, _d3_ace, + _d3_ade, + _d3_bcd, _d3_bce, + _d3_bde, _d3_cde; + mp::expansion< 96> _d4abcd, _d4abce, + _d4abde, _d4acde, + _d4bcde; + mp::expansion<960> _d5full; + + _OK = true; + + mp::expansion< 1 > _pa_zz_(_pa[ 2]); + mp::expansion< 1 > _pb_zz_(_pb[ 2]); + mp::expansion< 1 > _pc_zz_(_pc[ 2]); + mp::expansion< 1 > _pd_zz_(_pd[ 2]); + mp::expansion< 1 > _pe_zz_(_pe[ 2]); + + mp::expansion< 1 > _pa_qq_(_pa[ 3]); + mp::expansion< 1 > _pb_qq_(_pb[ 3]); + mp::expansion< 1 > _pc_qq_(_pc[ 3]); + mp::expansion< 1 > _pd_qq_(_pd[ 3]); + mp::expansion< 1 > _pe_qq_(_pe[ 3]); + + /*-------------------------------------- 2 x 2 minors */ + compute_det_2x2(_pa[ 0], _pa[ 1], + _pb[ 0], _pb[ 1], + _d2_ab_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pc[ 0], _pc[ 1], + _d2_ac_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pd[ 0], _pd[ 1], + _d2_ad_ ) ; + + compute_det_2x2(_pa[ 0], _pa[ 1], + _pe[ 0], _pe[ 1], + _d2_ae_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pc[ 0], _pc[ 1], + _d2_bc_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pd[ 0], _pd[ 1], + _d2_bd_ ) ; + + compute_det_2x2(_pb[ 0], _pb[ 1], + _pe[ 0], _pe[ 1], + _d2_be_ ) ; + + compute_det_2x2(_pc[ 0], _pc[ 1], + _pd[ 0], _pd[ 1], + _d2_cd_ ) ; + + compute_det_2x2(_pc[ 0], _pc[ 1], + _pe[ 0], _pe[ 1], + _d2_ce_ ) ; + + compute_det_2x2(_pd[ 0], _pd[ 1], + _pe[ 0], _pe[ 1], + _d2_de_ ) ; + + /*-------------------------------------- 3 x 3 minors */ + compute_det_3x3(_d2_bc_, _pa_zz_, + _d2_ac_, _pb_zz_, + _d2_ab_, _pc_zz_, + _d3_abc, +3) ; + + compute_det_3x3(_d2_bd_, _pa_zz_, + _d2_ad_, _pb_zz_, + _d2_ab_, _pd_zz_, + _d3_abd, +3) ; + + compute_det_3x3(_d2_be_, _pa_zz_, + _d2_ae_, _pb_zz_, + _d2_ab_, _pe_zz_, + _d3_abe, +3) ; + + compute_det_3x3(_d2_cd_, _pa_zz_, + _d2_ad_, _pc_zz_, + _d2_ac_, _pd_zz_, + _d3_acd, +3) ; + + compute_det_3x3(_d2_ce_, _pa_zz_, + _d2_ae_, _pc_zz_, + _d2_ac_, _pe_zz_, + _d3_ace, +3) ; + + compute_det_3x3(_d2_de_, _pa_zz_, + _d2_ae_, _pd_zz_, + _d2_ad_, _pe_zz_, + _d3_ade, +3) ; + + compute_det_3x3(_d2_cd_, _pb_zz_, + _d2_bd_, _pc_zz_, + _d2_bc_, _pd_zz_, + _d3_bcd, +3) ; + + compute_det_3x3(_d2_ce_, _pb_zz_, + _d2_be_, _pc_zz_, + _d2_bc_, _pe_zz_, + _d3_bce, +3) ; + + compute_det_3x3(_d2_de_, _pb_zz_, + _d2_be_, _pd_zz_, + _d2_bd_, _pe_zz_, + _d3_bde, +3) ; + + compute_det_3x3(_d2_de_, _pc_zz_, + _d2_ce_, _pd_zz_, + _d2_cd_, _pe_zz_, + _d3_cde, +3) ; + + /*-------------------------------------- 4 x 4 minors */ + unitary_det_4x4(_d3_cde, _d3_bde, + _d3_bce, _d3_bcd, + _d4bcde, +4) ; + + unitary_det_4x4(_d3_cde, _d3_ade, + _d3_ace, _d3_acd, + _d4acde, +4) ; + + unitary_det_4x4(_d3_bde, _d3_ade, + _d3_abe, _d3_abd, + _d4abde, +4) ; + + unitary_det_4x4(_d3_bce, _d3_ace, + _d3_abe, _d3_abc, + _d4abce, +4) ; + + unitary_det_4x4(_d3_bcd, _d3_acd, + _d3_abd, _d3_abc, + _d4abcd, +4) ; + + /*-------------------------------------- 5 x 5 result */ + compute_det_5x5(_d4bcde, _pa_qq_, + _d4acde, _pb_qq_, + _d4abde, _pc_qq_, + _d4abce, _pd_qq_, + _d4abcd, _pe_qq_, + _d5full, +4) ; + + /*-------------------------------------- leading det. */ + return mp::expansion_est(_d5full) ; + } + + __normal_call REAL_TYPE orient4d_i ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + __const_ptr(REAL_TYPE) _pd , + __const_ptr(REAL_TYPE) _pe , + bool_type &_OK + ) + { + /*--------------- orient4d predicate, "bound" version */ + ia_flt _aex, _aey, _aez , + _aeq, + _bex, _bey, _bez , + _beq, + _cex, _cey, _cez , + _ceq, + _dex, _dey, _dez , + _deq; + ia_flt _aexbey, _bexaey , + _aexcey, _cexaey , + _bexcey, _cexbey , + _cexdey, _dexcey , + _dexaey, _aexdey , + _bexdey, _dexbey ; + ia_flt _ab_, _bc_, _cd_, _da_, + _ac_, _bd_; + ia_flt _abc, _bcd, _cda, _dab; + ia_flt _sgn; + + ia_rnd _rnd; // up rounding! + + _aex.from_sub(_pa[0], _pe[0]) ; // coord. diff. + _aey.from_sub(_pa[1], _pe[1]) ; + _aez.from_sub(_pa[2], _pe[2]) ; + _aeq.from_sub(_pa[3], _pe[3]) ; + + _bex.from_sub(_pb[0], _pe[0]) ; + _bey.from_sub(_pb[1], _pe[1]) ; + _bez.from_sub(_pb[2], _pe[2]) ; + _beq.from_sub(_pb[3], _pe[3]) ; + + _cex.from_sub(_pc[0], _pe[0]) ; + _cey.from_sub(_pc[1], _pe[1]) ; + _cez.from_sub(_pc[2], _pe[2]) ; + _ceq.from_sub(_pc[3], _pe[3]) ; + + _dex.from_sub(_pd[0], _pe[0]) ; + _dey.from_sub(_pd[1], _pe[1]) ; + _dez.from_sub(_pd[2], _pe[2]) ; + _deq.from_sub(_pd[3], _pe[3]) ; + + _aexbey = _aex * _bey ; // 2 x 2 minors + _bexaey = _bex * _aey ; + _ab_ = _aexbey - _bexaey ; + + _bexcey = _bex * _cey; + _cexbey = _cex * _bey; + _bc_ = _bexcey - _cexbey ; + + _cexdey = _cex * _dey; + _dexcey = _dex * _cey; + _cd_ = _cexdey - _dexcey ; + + _dexaey = _dex * _aey; + _aexdey = _aex * _dey; + _da_ = _dexaey - _aexdey ; + + _aexcey = _aex * _cey; + _cexaey = _cex * _aey; + _ac_ = _aexcey - _cexaey ; + + _bexdey = _bex * _dey; + _dexbey = _dex * _bey; + _bd_ = _bexdey - _dexbey ; + + _abc = // 3 x 3 minors + _aez * _bc_ - _bez * _ac_ + + _cez * _ab_ ; + + _bcd = + _bez * _cd_ - _cez * _bd_ + + _dez * _bc_ ; + + _cda = + _cez * _da_ + _dez * _ac_ + + _aez * _cd_ ; + + _dab = + _dez * _ab_ + _aez * _bd_ + + _bez * _da_ ; + + _sgn = // 4 x 4 result + _deq * _abc - _ceq * _dab + + _beq * _cda - _aeq * _bcd ; + + _OK = + _sgn.lo() >= (REAL_TYPE)0. + ||_sgn.up() <= (REAL_TYPE)0.; + + return ( _sgn.mid() ) ; + } + + __normal_call REAL_TYPE orient4d_f ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + __const_ptr(REAL_TYPE) _pd , + __const_ptr(REAL_TYPE) _pe , + bool_type &_OK + ) + { + /*--------------- orient4d predicate, "float" version */ + REAL_TYPE static const _ER = + + 13. * std::pow(mp::_epsilon, 1) ; + + REAL_TYPE _aex, _aey, _aez , + _aeq, + _bex, _bey, _bez , + _beq, + _cex, _cey, _cez , + _ceq, + _dex, _dey, _dez , + _deq; + REAL_TYPE _aexbey, _bexaey , + _aexcey, _cexaey , + _bexcey, _cexbey , + _cexdey, _dexcey , + _dexaey, _aexdey , + _bexdey, _dexbey ; + REAL_TYPE _ab_, _bc_, _cd_, _da_, + _ac_, _bd_; + REAL_TYPE _abc, _bcd, _cda, _dab; + + REAL_TYPE _AEZ, _BEZ, _CEZ, _DEZ; + REAL_TYPE _AEQ, _BEQ, _CEQ, _DEQ; + REAL_TYPE _AEXBEY, _BEXAEY , + _CEXAEY, _AEXCEY , + _BEXCEY, _CEXBEY , + _CEXDEY, _DEXCEY , + _DEXAEY, _AEXDEY , + _BEXDEY, _DEXBEY ; + REAL_TYPE _AB_, _BC_, _CD_, _DA_, + _AC_, _BD_; + REAL_TYPE _ABC, _BCD, _CDA, _DAB; + + REAL_TYPE _sgn, _FT ; + + _aex = _pa [0] - _pe [0] ; // coord. diff. + _aey = _pa [1] - _pe [1] ; + _aez = _pa [2] - _pe [2] ; + _aeq = _pa [3] - _pe [3] ; + + _AEZ = std::abs (_aez) ; + _AEQ = std::abs (_aeq) ; + + _bex = _pb [0] - _pe [0] ; + _bey = _pb [1] - _pe [1] ; + _bez = _pb [2] - _pe [2] ; + _beq = _pb [3] - _pe [3] ; + + _BEZ = std::abs (_bez) ; + _BEQ = std::abs (_beq) ; + + _cex = _pc [0] - _pe [0] ; + _cey = _pc [1] - _pe [1] ; + _cez = _pc [2] - _pe [2] ; + _ceq = _pc [3] - _pe [3] ; + + _CEZ = std::abs (_cez) ; + _CEQ = std::abs (_ceq) ; + + _dex = _pd [0] - _pe [0] ; + _dey = _pd [1] - _pe [1] ; + _dez = _pd [2] - _pe [2] ; + _deq = _pd [3] - _pe [3] ; + + _DEZ = std::abs (_dez) ; + _DEQ = std::abs (_deq) ; + + _aexbey = _aex * _bey; // 2 x 2 minors + _bexaey = _bex * _aey; + _ab_ = _aexbey - _bexaey ; + + _AEXBEY = std::abs (_aexbey) ; + _BEXAEY = std::abs (_bexaey) ; + _AB_ = _AEXBEY + _BEXAEY ; + + _bexcey = _bex * _cey; + _cexbey = _cex * _bey; + _bc_ = _bexcey - _cexbey ; + + _BEXCEY = std::abs (_bexcey) ; + _CEXBEY = std::abs (_cexbey) ; + _BC_ = _BEXCEY + _CEXBEY ; + + _cexdey = _cex * _dey; + _dexcey = _dex * _cey; + _cd_ = _cexdey - _dexcey ; + + _CEXDEY = std::abs (_cexdey) ; + _DEXCEY = std::abs (_dexcey) ; + _CD_ = _CEXDEY + _DEXCEY ; + + _dexaey = _dex * _aey; + _aexdey = _aex * _dey; + _da_ = _dexaey - _aexdey ; + + _DEXAEY = std::abs (_dexaey) ; + _AEXDEY = std::abs (_aexdey) ; + _DA_ = _DEXAEY + _AEXDEY ; + + _aexcey = _aex * _cey; + _cexaey = _cex * _aey; + _ac_ = _aexcey - _cexaey ; + + _AEXCEY = std::abs (_aexcey) ; + _CEXAEY = std::abs (_cexaey) ; + _AC_ = _AEXCEY + _CEXAEY ; + + _bexdey = _bex * _dey; + _dexbey = _dex * _bey; + _bd_ = _bexdey - _dexbey ; + + _BEXDEY = std::abs (_bexdey) ; + _DEXBEY = std::abs (_dexbey) ; + _BD_ = _BEXDEY + _DEXBEY ; + + _abc = // 3 x 3 minors + _aez * _bc_ - _bez * _ac_ + + _cez * _ab_ ; + _ABC = + _AEZ * _BC_ + _BEZ * _AC_ + + _CEZ * _AB_ ; + + _bcd = + _bez * _cd_ - _cez * _bd_ + + _dez * _bc_ ; + _BCD = + _BEZ * _CD_ + _CEZ * _BD_ + + _DEZ * _BC_ ; + + _cda = + _cez * _da_ + _dez * _ac_ + + _aez * _cd_ ; + _CDA = + _CEZ * _DA_ + _DEZ * _AC_ + + _AEZ * _CD_ ; + + _dab = + _dez * _ab_ + _aez * _bd_ + + _bez * _da_ ; + _DAB = + _DEZ * _AB_ + _AEZ * _BD_ + + _BEZ * _DA_ ; + + _FT = // roundoff tol + _DEQ * _ABC + _CEQ * _DAB + + _BEQ * _CDA + _AEQ * _BCD ; + + _FT *= _ER ; + + _sgn = // 4 x 4 result + _deq * _abc - _ceq * _dab + + _beq * _cda - _aeq * _bcd ; + + _OK = + _sgn > _FT || _sgn < -_FT ; + + return ( _sgn ) ; + } diff --git a/predicate/predicate_k.hpp b/predicate/predicate_k.hpp index d28493f..5ab6be0 100644 --- a/predicate/predicate_k.hpp +++ b/predicate/predicate_k.hpp @@ -47,20 +47,24 @@ * how they can obtain it for free, then you are not * required to make any arrangement with me.) * - * Disclaimer: Neither I nor: Columbia University, The - * Massachusetts Institute of Technology, The - * University of Sydney, nor The National Aeronautics - * and Space Administration warrant this code in any - * way whatsoever. This code is provided "as-is" to be - * used at your own risk. + * Disclaimer: Neither I nor THE CONTRIBUTORS warrant + * this code in any way whatsoever. This code is + * provided "as-is" to be used at your own risk. + * + * THE CONTRIBUTORS include: + * (a) The University of Sydney + * (b) The Massachusetts Institute of Technology + * (c) Columbia University + * (d) The National Aeronautics & Space Administration + * (e) Los Alamos National Laboratory * -------------------------------------------------------- * - * Last updated: 15 April, 2020 + * Last updated: 11 May., 2024 * * Copyright 2020-- * Darren Engwirda - * de2363@columbia.edu + * d.engwirda@gmail.com * https://github.com/dengwirda/ * -------------------------------------------------------- @@ -72,7 +76,7 @@ # define __PREDICATE_K__ # define USE_KERNEL_FLTPOINT -// define USE_KERNEL_INTERVAL +# define USE_KERNEL_INTERVAL namespace geompred { @@ -84,14 +88,19 @@ enum _kernel { ORIENT2D_f, ORIENT2D_i, ORIENT2D_e , ORIENT3D_f, ORIENT3D_i, ORIENT3D_e , + ORIENT4D_f, ORIENT4D_i, ORIENT4D_e , BISECT2D_f, BISECT2D_i, BISECT2D_e , BISECT2W_f, BISECT2W_i, BISECT2W_e , BISECT3D_f, BISECT3D_i, BISECT3D_e , BISECT3W_f, BISECT3W_i, BISECT3W_e , + BISECT4D_f, BISECT4D_i, BISECT4D_e , + BISECT4W_f, BISECT4W_i, BISECT4W_e , INBALL2D_f, INBALL2D_i, INBALL2D_e , INBALL2W_f, INBALL2W_i, INBALL2W_e , INBALL3D_f, INBALL3D_i, INBALL3D_e , INBALL3W_f, INBALL3W_i, INBALL3W_e , + INBALL4D_f, INBALL4D_i, INBALL4D_e , + INBALL4W_f, INBALL4W_i, INBALL4W_e , LASTKERNEL } ; size_t _nn_calls[LASTKERNEL] = {0} ; @@ -186,6 +195,50 @@ return (REAL_TYPE) +0.0E+00; } + __inline_call REAL_TYPE orient4d ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc , + __const_ptr(REAL_TYPE) _pd , + __const_ptr(REAL_TYPE) _pe + ) + { + /*------------ orient4d predicate, "filtered" version */ + REAL_TYPE _rr; + bool_type _OK; + + # ifdef USE_KERNEL_FLTPOINT + _nn_calls[ORIENT4D_f] += +1; + + _rr = orient4d_f( // "float" kernel + _pa, _pb, _pc, _pd, _pe, _OK + ) ; + + if (_OK && std::isnormal(_rr)) + return _rr ; + # endif + + # ifdef USE_KERNEL_INTERVAL + _nn_calls[ORIENT4D_i] += +1; + + _rr = orient4d_i( // "bound" kernel + _pa, _pb, _pc, _pd, _pe, _OK + ) ; + + if (_OK) return _rr ; + # endif + + _nn_calls[ORIENT4D_e] += +1; + + _rr = orient4d_e( // "exact" kernel + _pa, _pb, _pc, _pd, _pe, _OK + ) ; + + if (_OK) return _rr ; + + return (REAL_TYPE) +0.0E+00; + } + __inline_call REAL_TYPE bisect2d ( __const_ptr(REAL_TYPE) _pa , __const_ptr(REAL_TYPE) _pb , @@ -368,6 +421,97 @@ } } + __inline_call REAL_TYPE bisect4d ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc + ) + { + /*------------ bisect4d predicate, "filtered" version */ + REAL_TYPE _rr; + bool_type _OK; + + # ifdef USE_KERNEL_FLTPOINT + _nn_calls[BISECT4D_f] += +1; + + _rr = bisect4d_f( // "float" kernel + _pa, _pb, _pc, _OK + ) ; + + if (_OK && std::isnormal(_rr)) + return _rr ; + # endif + + # ifdef USE_KERNEL_INTERVAL + _nn_calls[BISECT4D_i] += +1; + + _rr = bisect4d_i( // "bound" kernel + _pa, _pb, _pc, _OK + ) ; + + if (_OK) return _rr ; + # endif + + _nn_calls[BISECT4D_e] += +1; + + _rr = bisect4d_e( // "exact" kernel + _pa, _pb, _pc, _OK + ) ; + + if (_OK) return _rr ; + + return (REAL_TYPE) +0.0E+00; + } + + __inline_call REAL_TYPE bisect4w ( + __const_ptr(REAL_TYPE) _pa , + __const_ptr(REAL_TYPE) _pb , + __const_ptr(REAL_TYPE) _pc + ) + { + /*------------ bisect3w predicate, "filtered" version */ + if (_pa [ 4] == _pb [ 4] ) + { // equal weights, do bisect3d + return bisect4d(_pa, _pb, _pc) ; + } + else + { + REAL_TYPE _rr; // given weights, full kernel + bool_type _OK; + + # ifdef USE_KERNEL_FLTPOINT + _nn_calls[BISECT4W_f] += +1; + + _rr = bisect4w_f( // "float" kernel + _pa, _pb, _pc, _OK + ) ; + + if (_OK && std::isnormal(_rr)) + return _rr ; + # endif + + # ifdef USE_KERNEL_INTERVAL + _nn_calls[BISECT4W_i] += +1; + + _rr = bisect4w_i( // "bound" kernel + _pa, _pb, _pc, _OK + ) ; + + if (_OK) return _rr ; + # endif + + _nn_calls[BISECT4W_e] += +1; + + _rr = bisect4w_e( // "exact" kernel + _pa, _pb, _pc, _OK + ) ; + + if (_OK) return _rr ; + + return (REAL_TYPE) +0.0E+00; + } + } + __inline_call REAL_TYPE inball2d ( __const_ptr(REAL_TYPE) _pa , __const_ptr(REAL_TYPE) _pb ,