MRPrecisePredicates3Poly.h

Go to the documentation of this file.

#pragma once
 
#include "MRPrecisePredicates3.h"
#include "MRInt64Mul128.h"
#include "MRMacros.h"
#include "MRSparsePolynomial.h"
#include <algorithm>
#include <array>
#include <cassert>
#include <limits>
#include <type_traits>
#include <vector>
 
namespace MR
{


namespace detail
{


[[nodiscard]] MR_FORCE_INLINE FastInt128 orient3dVolume( const Vector3i & a, const Vector3i & b, const Vector3i & c, const Vector3i & d )
{
    const Vector3i64 x( a - d );
    const Vector3i64 y( b - d );
    const Vector3i64 z( c - d );
 
    return
        Int64Mul128( x.x ) * Int64Mul128( y.y * z.z - y.z * z.y )
     -  Int64Mul128( x.y ) * Int64Mul128( y.x * z.z - y.z * z.x )
     +  Int64Mul128( x.z ) * Int64Mul128( y.x * z.y - y.y * z.x );
}

struct PointDegree
{
    Vector3i pt;
    std::int64_t d = 0; 

    [[nodiscard]] std::int64_t degree( int i ) const
    {
        constexpr std::int64_t w[3] = { 9, 3, 1 };
        return w[i] * d;
    }
};

constexpr int cMaxPerturbedPoints = 16;

template <int N = cMaxPerturbedPoints, bool Bounded = false>
MR_FORCE_INLINE int getPointDegrees( const PreciseVertCoords * vs, const std::int64_t * radices, int n, PointDegree * res,
    std::int64_t * places = nullptr )
{
    struct VertN
    {
        VertId v;
        int n = 0;
    };
    assert( n > 0 && n <= N );
    std::array<VertN, N> as;
    for ( int i = 0; i < n; ++i )
        as[i] = { vs[i].id, i };
    std::sort( as.begin(), as.begin() + n, []( const auto & a, const auto & b ) { return a.v < b.v; } );
 
    std::int64_t d = 1;
    int numIds = 1;
    if ( places )
        places[0] = d;
    for ( int i = 0; i < n; ++i )
    {
        const auto k = as[i].n;
        if constexpr ( Bounded )
        {
            if ( d > std::numeric_limits<std::int64_t>::max() / radices[k] )
            {
                for ( int j = i; j < n; ++j )
                    res[as[j].n] = { vs[as[j].n].pt, 0 };
                return numIds - 1;
            }
        }
        res[k] = { vs[k].pt, d };
        if ( i + 1 == n )
            break;
        const auto r = radices[k];
        if ( as[i].v < as[i+1].v ) 
        {
            assert( d <= std::numeric_limits<std::int64_t>::max() / r );
            d *= r;
            if ( places )
                places[numIds] = d;
            ++numIds;
        }
        else
            assert( r == radices[as[i+1].n] );
    }
    return numIds;
}

constexpr std::int64_t cMaxPolyDAll = std::numeric_limits<std::int64_t>::max();

constexpr std::int64_t cMaxPolyDTriTri   = 15'943'959; 
constexpr std::int64_t cMaxPolyDTriPlane = 430'486'893; 

template <int N, bool Bounded>
class PerturbationStagesOf
{
public:
    MR_FORCE_INLINE PerturbationStagesOf( const PreciseVertCoords * vs, const std::int64_t * radices, int n )
    {
        assert( n <= N );
        numStages_ = getPointDegrees<N, Bounded>( vs, radices, n, points_.data(), places_.data() );
    }

    [[nodiscard]] int numStages() const { return numStages_; }

    [[nodiscard]] std::int64_t place( int s ) const { return places_[s - 1]; }

    [[nodiscard]] const PointDegree & operator[]( int i ) const { return points_[i]; }

    [[nodiscard]] PointDegree at( int i, int s ) const { return points_[i].d <= place( s ) ? points_[i] : PointDegree{ points_[i].pt, 0 }; }

    [[nodiscard]] bool startsAt( int i, int s ) const { return points_[i].d == place( s ); }

    template <typename F>
    [[nodiscard]] auto firstStage( F && f, std::int64_t maxDeg = cMaxPolyDAll ) const -> std::invoke_result_t<F &, int>
    {
        for ( int s = 1; s <= numStages_ && place( s ) <= maxDeg; ++s )
            if ( auto res = f( s ) )
                return res;
        return {};
    }
 
private:
    std::array<PointDegree, N> points_;
    std::array<std::int64_t, N> places_;
    int numStages_ = 0;
};

using PerturbationStages = PerturbationStagesOf<cMaxPerturbedPoints, false>;

template <std::int64_t M>
using VolumePoly = SparsePolynomial<FastInt128, std::int64_t, M>;

template <std::int64_t M>
[[nodiscard]] VolumePoly<M> orient3dPoly( const PointDegree & a, const PointDegree & b, const PointDegree & c, const PointDegree & d )
{
    struct Element
    {
        typename VolumePoly<M>::Term terms[3];
        int numTerms = 0;
    };
    Element m[3][3];
    const PointDegree * rows[3] = { &a, &b, &c };
    for ( int r = 0; r < 3; ++r )
        for ( int col = 0; col < 3; ++col )
        {
            auto & e = m[r][col];
            if ( const std::int64_t delta = std::int64_t( rows[r]->pt[col] ) - d.pt[col] )
                e.terms[e.numTerms++] = { 0, delta };
            if ( rows[r]->d )
                e.terms[e.numTerms++] = { rows[r]->degree( col ), 1 };
            if ( d.d )
                e.terms[e.numTerms++] = { d.degree( col ), -1 };
        }

    constexpr int perms[6][4] = { { 0, 1, 2, 1 }, { 1, 2, 0, 1 }, { 2, 0, 1, 1 }, { 0, 2, 1, -1 }, { 2, 1, 0, -1 }, { 1, 0, 2, -1 } };
    std::vector<typename VolumePoly<M>::Term> terms;
    terms.reserve( 6 * 26 );
    for ( const auto & p : perms )
    {
        const Element & e0 = m[0][p[0]];
        const Element & e1 = m[1][p[1]];
        const Element & e2 = m[2][p[2]];
        for ( int k0 = 0; k0 < e0.numTerms; ++k0 )
            for ( int k1 = 0; k1 < e1.numTerms; ++k1 )
            {
                const std::int64_t c01 = p[3] * std::int64_t( e0.terms[k0].second ) * std::int64_t( e1.terms[k1].second );
                const std::int64_t d01 = e0.terms[k0].first + e1.terms[k1].first;
                const int last2 = d.d && k0 + 1 == e0.numTerms && k1 + 1 == e1.numTerms ? e2.numTerms - 1 : e2.numTerms;
                for ( int k2 = 0; k2 < last2; ++k2 )
                    if ( const auto deg = d01 + e2.terms[k2].first; deg <= M )
                        terms.emplace_back( deg, Int64Mul128( c01 ) * Int64Mul128( std::int64_t( e2.terms[k2].second ) ) );
            }
    }
    return VolumePoly<M>::fromUnsortedTerms( std::move( terms ) );
}
 
}

 
}