Skip to content

Commit af729f3

Browse files
FedrFedr
authored
SparsePolynomial: limit maximum degree to store inside (#4795)
1 parent 8079fd4 commit af729f3

File tree

3 files changed

+74
-42
lines changed

3 files changed

+74
-42
lines changed

‎source/MRMesh/MRPrecisePredicates2.cpp

Lines changed: 13 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -39,14 +39,13 @@ std::array<PointDegree, 6> getPointDegrees( const std::array<PreciseVertCoords2,
3939
return res;
4040
}
4141

42-
SparsePolynomial<FastInt128> ccwPoly(
43-
const PointDegree & a,
44-
const PointDegree & b,
45-
const PointDegree & c,
42+
// 840 was found experimentally for segmentIntersectionOrder with all 6 points have equal coordinates (but different ids);
43+
// if it is not enough then we will get assert violation inside poly.isPositive(), and increase the value
44+
using Poly = SparsePolynomial<FastInt128, int, 840>;
45+
46+
Poly ccwPoly( const PointDegree & a, const PointDegree & b, const PointDegree & c,
4647
int db ) // degree.x = degree.y * db
4748
{
48-
using Poly = SparsePolynomial<FastInt128>;
49-
5049
const Poly xx( a.pt.x - c.pt.x, a.d * db, 1, c.d * db, -1 );
5150
const Poly xy( a.pt.y - c.pt.y, a.d , 1, c.d , -1 );
5251
const Poly yx( b.pt.x - c.pt.x, b.d * db, 1, c.d * db, -1 );
@@ -56,16 +55,6 @@ SparsePolynomial<FastInt128> ccwPoly(
5655
return det;
5756
}
5857

59-
bool isPositive( const SparsePolynomial<FastInt128>& poly )
60-
{
61-
const auto & mapDegToCf = poly.get();
62-
if ( !mapDegToCf.empty() )
63-
return mapDegToCf.begin()->second > 0;
64-
65-
assert (false);
66-
return false;
67-
}
68-
6958
} // anonymous namespace
7059

7160
// see https://arxiv.org/pdf/math/9410209 Table 4-i:
@@ -302,17 +291,21 @@ bool segmentIntersectionOrder( const std::array<PreciseVertCoords2, 6> & vs )
302291
// ( ccw(sa,s[0])-ccw(sa,ds[1]) ) * ( ccw(sb,s[0])-ccw(sb,ds[1]) )
303292
const auto polySaOrg = ccwPoly( ds[2], ds[3], ds[0], 3 );
304293
const auto polySaDest = ccwPoly( ds[2], ds[3], ds[1], 3 );
305-
assert( isPositive( polySaOrg ) != isPositive( polySaDest ) );
294+
assert( !polySaOrg.empty() || !polySaDest.empty() );
295+
assert( polySaOrg.empty() || polySaDest.empty() || polySaOrg.isPositive() != polySaDest.isPositive() );
296+
const bool posSaOrg = polySaOrg.empty() ? !polySaDest.isPositive() : polySaOrg.isPositive();
306297

307298
const auto polySbOrg = ccwPoly( ds[4], ds[5], ds[0], 3 );
308299
const auto polySbDest = ccwPoly( ds[4], ds[5], ds[1], 3 );
309-
assert( isPositive( polySbOrg ) != isPositive( polySbDest ) );
300+
assert( !polySbOrg.empty() || !polySbDest.empty() );
301+
assert( polySbOrg.empty() || polySbDest.empty() || polySbOrg.isPositive() != polySbDest.isPositive() );
302+
const bool posSbOrg = polySbOrg.empty() ? !polySbDest.isPositive() : polySbOrg.isPositive();
310303

311304
auto nom = polySaOrg * polySbDest;
312305
nom -= polySbOrg * polySaDest;
313306

314-
bool res = isPositive( nom );
315-
if ( isPositive( polySaOrg ) != isPositive( polySbOrg ) ) // denominator is negative
307+
bool res = nom.isPositive();
308+
if ( posSaOrg != posSbOrg ) // denominator is negative
316309
res = !res;
317310
return res;
318311
}

‎source/MRMesh/MRSparsePolynomial.h

Lines changed: 46 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -5,17 +5,19 @@
55
namespace MR
66
{
77

8-
template <typename C, typename D>
8+
template <typename C, typename D, D M>
99
class SparsePolynomial;
10-
template <typename C, typename D>
11-
SparsePolynomial<C,D> operator *( const SparsePolynomial<C,D>& a, const SparsePolynomial<C,D>& b );
10+
template <typename C, typename D, D M>
11+
SparsePolynomial<C,D,M> operator *( const SparsePolynomial<C,D,M>& a, const SparsePolynomial<C,D,M>& b );
1212

1313
/// The class to store a polynomial with a large number of zero coefficient (only non-zeros are stored in std::map)
1414
/// \tparam C - type of coefficients
1515
/// \tparam D - type of degrees
16-
template <typename C, typename D = int>
16+
/// \tparam M - maximum degree to store in the polynomial
17+
template <typename C, typename D, D M>
1718
class SparsePolynomial
1819
{
20+
static_assert( M > 0 );
1921
public:
2022
/// constructs zero polynomial
2123
SparsePolynomial() = default;
@@ -29,8 +31,14 @@ class SparsePolynomial
2931
/// constructs polynomial c0 + c1*x^d1 + c2*x^d2
3032
SparsePolynomial( C c0, D d1, C c1, D d2, C c2 );
3133

34+
/// returns true if no single polynomial coefficient is defined
35+
[[nodiscard]] bool empty() const { return map_.empty(); }
36+
37+
/// returns true if the coefficient for the smallest not-zero degress is positive
38+
[[nodiscard]] bool isPositive() const;
39+
3240
/// gets read-only access to all not-zero coefficients
33-
const std::map<D, C> & get() const { return map_; }
41+
[[nodiscard]] const std::map<D, C> & get() const { return map_; }
3442

3543
SparsePolynomial& operator +=( const SparsePolynomial& b );
3644
SparsePolynomial& operator -=( const SparsePolynomial& b );
@@ -40,27 +48,31 @@ class SparsePolynomial
4048
std::map<D, C> map_; // degree -> not-zero coefficient
4149
};
4250

43-
template <typename C, typename D>
44-
SparsePolynomial<C,D>::SparsePolynomial( std::map<D, C> && m ) : map_( std::move( m ) )
51+
template <typename C, typename D, D M>
52+
SparsePolynomial<C,D,M>::SparsePolynomial( std::map<D, C> && m ) : map_( std::move( m ) )
4553
{
4654
#ifndef NDEBUG
4755
for ( const auto & [deg, cf] : map_ )
56+
{
57+
assert( deg <= M );
4858
assert( cf != 0 );
59+
}
4960
#endif
5061
}
5162

52-
template <typename C, typename D>
53-
SparsePolynomial<C,D>::SparsePolynomial( C c0, D d1, C c1 )
63+
template <typename C, typename D, D M>
64+
SparsePolynomial<C,D,M>::SparsePolynomial( C c0, D d1, C c1 )
5465
{
5566
assert( c1 != 0 );
5667
assert( d1 != 0 );
5768
if ( c0 != 0 )
5869
map_[D(0)] = c0;
59-
map_[d1] = c1;
70+
if ( d1 <= M )
71+
map_[d1] = c1;
6072
}
6173

62-
template <typename C, typename D>
63-
SparsePolynomial<C,D>::SparsePolynomial( C c0, D d1, C c1, D d2, C c2 )
74+
template <typename C, typename D, D M>
75+
SparsePolynomial<C,D,M>::SparsePolynomial( C c0, D d1, C c1, D d2, C c2 )
6476
{
6577
assert( c1 != 0 );
6678
assert( d1 != 0 );
@@ -69,12 +81,24 @@ SparsePolynomial<C,D>::SparsePolynomial( C c0, D d1, C c1, D d2, C c2 )
6981
assert( d1 != d2 );
7082
if ( c0 != 0 )
7183
map_[D(0)] = c0;
72-
map_[d1] = c1;
73-
map_[d2] = c2;
84+
if ( d1 <= M )
85+
map_[d1] = c1;
86+
if ( d2 <= M )
87+
map_[d2] = c2;
88+
}
89+
90+
template <typename C, typename D, D M>
91+
bool SparsePolynomial<C,D,M>::isPositive() const
92+
{
93+
if ( !map_.empty() )
94+
return map_.begin()->second > 0;
95+
96+
assert (false);
97+
return false;
7498
}
7599

76-
template <typename C, typename D>
77-
SparsePolynomial<C,D>& SparsePolynomial<C,D>::operator +=( const SparsePolynomial& b )
100+
template <typename C, typename D, D M>
101+
SparsePolynomial<C,D,M>& SparsePolynomial<C,D,M>::operator +=( const SparsePolynomial& b )
78102
{
79103
for ( const auto & [degB, cfB] : b.map_ )
80104
{
@@ -92,8 +116,8 @@ SparsePolynomial<C,D>& SparsePolynomial<C,D>::operator +=( const SparsePolynomia
92116
return * this;
93117
}
94118

95-
template <typename C, typename D>
96-
SparsePolynomial<C,D>& SparsePolynomial<C,D>::operator -=( const SparsePolynomial& b )
119+
template <typename C, typename D, D M>
120+
SparsePolynomial<C,D,M>& SparsePolynomial<C,D,M>::operator -=( const SparsePolynomial& b )
97121
{
98122
for ( const auto & [degB, cfB] : b.map_ )
99123
{
@@ -111,8 +135,8 @@ SparsePolynomial<C,D>& SparsePolynomial<C,D>::operator -=( const SparsePolynomia
111135
return * this;
112136
}
113137

114-
template <typename C, typename D>
115-
SparsePolynomial<C,D> operator *( const SparsePolynomial<C,D>& a, const SparsePolynomial<C,D>& b )
138+
template <typename C, typename D, D M>
139+
[[nodiscard]] SparsePolynomial<C,D,M> operator *( const SparsePolynomial<C,D,M>& a, const SparsePolynomial<C,D,M>& b )
116140
{
117141
std::map<D,C> resMap;
118142
for ( const auto & [degA, cfA] : a.map_ )
@@ -122,6 +146,8 @@ SparsePolynomial<C,D> operator *( const SparsePolynomial<C,D>& a, const SparsePo
122146
{
123147
assert( cfB != 0 );
124148
const auto deg = degA + degB;
149+
if ( deg > M )
150+
break;
125151
const auto cf = cfA * cfB;
126152
auto [it,inserted] = resMap.insert( { deg, cf } );
127153
if ( !inserted )

‎source/MRTest/MRPrecisePredicatesTests.cpp

Lines changed: 15 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -121,17 +121,30 @@ TEST( MRMesh, PrecisePredicates2More )
121121

122122
TEST( MRMesh, PrecisePredicates2FullDegen )
123123
{
124-
std::array<PreciseVertCoords2, 4> vs =
124+
std::array<PreciseVertCoords2, 6> vs =
125125
{
126126
PreciseVertCoords2{ 0_v, Vector2i( 0, 0 ) },
127127
PreciseVertCoords2{ 1_v, Vector2i( 0, 0 ) },
128128
PreciseVertCoords2{ 2_v, Vector2i( 0, 0 ) },
129-
PreciseVertCoords2{ 3_v, Vector2i( 0, 0 ) }
129+
PreciseVertCoords2{ 3_v, Vector2i( 0, 0 ) },
130+
PreciseVertCoords2{ 4_v, Vector2i( 0, 0 ) },
131+
PreciseVertCoords2{ 5_v, Vector2i( 0, 0 ) }
130132
};
131133

132134
EXPECT_FALSE( doSegmentSegmentIntersect( { vs[0], vs[1], vs[2], vs[3] } ).doIntersect );
133135
EXPECT_TRUE( doSegmentSegmentIntersect( { vs[0], vs[2], vs[1], vs[3] } ).doIntersect );
134136
EXPECT_FALSE( doSegmentSegmentIntersect( { vs[0], vs[3], vs[1], vs[2] } ).doIntersect );
137+
138+
// test that maximum degree in segmentIntersectionOrder can cope with most degenerate situation possible
139+
do
140+
{
141+
if ( doSegmentSegmentIntersect( { vs[0], vs[1], vs[2], vs[3] } ).doIntersect
142+
&& doSegmentSegmentIntersect( { vs[0], vs[1], vs[4], vs[5] } ).doIntersect )
143+
{
144+
(void)segmentIntersectionOrder( { vs[0], vs[1], vs[2], vs[3], vs[4], vs[5] } );
145+
}
146+
}
147+
while ( std::next_permutation( vs.begin(), vs.end(), []( const auto & l, const auto & r ) { return l.id < r.id; } ) );
135148
}
136149

137150
TEST( MRMesh, PrecisePredicates2PartialDegen )

0 commit comments

Comments
 (0)