GCC Code Coverage Report

Directory: ./
File: lib/geogram/numerics/expansion_nt.cpp
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 36 99 36.4%
Functions: 7 17 41.2%
Branches: 11 14 78.6%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/numerics/expansion_nt.h>
41
42 namespace GEO {
43
44 ✗ expansion_nt& expansion_nt::operator+= (const expansion_nt& rhs) {
45 index_t e_capa = expansion::sum_capacity(rep(), rhs.rep());
46 expansion* e = expansion::new_expansion_on_heap(e_capa);
47 ✗ e->assign_sum(rep(), rhs.rep());
48 cleanup();
49 ✗ rep_ = e;
50 ✗ return *this;
51 }
52
53 ✗ expansion_nt& expansion_nt::operator+= (double rhs) {
54 index_t e_capa = expansion::sum_capacity(rep(), rhs);
55
56 // TODO: optimized in-place version to be used
57 // if(!shared() && e_capa < rep().capacity())
58
59 expansion* e = expansion::new_expansion_on_heap(e_capa);
60 ✗ e->assign_sum(rep(), rhs);
61 cleanup();
62 ✗ rep_ = e;
63 ✗ return *this;
64 }
65
66 ✗ expansion_nt& expansion_nt::operator-= (const expansion_nt& rhs) {
67 index_t e_capa = expansion::diff_capacity(rep(), rhs.rep());
68 expansion* e = expansion::new_expansion_on_heap(e_capa);
69 ✗ e->assign_diff(rep(), rhs.rep());
70 cleanup();
71 ✗ rep_ = e;
72 ✗ return *this;
73 }
74
75 ✗ expansion_nt& expansion_nt::operator-= (double rhs) {
76 index_t e_capa = expansion::diff_capacity(rep(), rhs);
77
78 // TODO: optimized in-place version to be used
79 // if(!shared() && e_capa < rep().capacity())
80
81 expansion* e = expansion::new_expansion_on_heap(e_capa);
82 ✗ e->assign_diff(rep(), rhs);
83 cleanup();
84 ✗ rep_ = e;
85 ✗ return *this;
86 }
87
88 ✗ expansion_nt& expansion_nt::operator*= (const expansion_nt& rhs) {
89 index_t e_capa = expansion::product_capacity(rep(), rhs.rep());
90 expansion* e = expansion::new_expansion_on_heap(e_capa);
91 ✗ e->assign_product(rep(), rhs.rep());
92 cleanup();
93 ✗ rep_ = e;
94 ✗ return *this;
95 }
96
97 ✗ expansion_nt& expansion_nt::operator*= (double rhs) {
98 index_t e_capa = expansion::product_capacity(rep(), rhs);
99
100 // TODO: optimized in-place version to be used
101 // if(!shared() && e_capa < rep().capacity())
102
103 expansion* e = expansion::new_expansion_on_heap(e_capa);
104 ✗ e->assign_product(rep(), rhs);
105 cleanup();
106 ✗ rep_ = e;
107 ✗ return *this;
108 }
109
110 /************************************************************************/
111
112 437550 expansion_nt expansion_nt::operator+ (const expansion_nt& rhs) const {
113 expansion* e = expansion::new_expansion_on_heap(
114 expansion::sum_capacity(rep(), rhs.rep())
115 );
116 437550 e->assign_sum(rep(), rhs.rep());
117 437550 return expansion_nt(e);
118 }
119
120 1482723 expansion_nt expansion_nt::operator- (const expansion_nt& rhs) const {
121 expansion* e = expansion::new_expansion_on_heap(
122 expansion::diff_capacity(rep(), rhs.rep())
123 );
124 1482723 e->assign_diff(rep(), rhs.rep());
125 1482723 return expansion_nt(e);
126 }
127
128 1517093 expansion_nt expansion_nt::operator* (const expansion_nt& rhs) const {
129 expansion* e = expansion::new_expansion_on_heap(
130 expansion::product_capacity(rep(), rhs.rep())
131 );
132 1517093 e->assign_product(rep(), rhs.rep());
133 1517093 return expansion_nt(e);
134 }
135
136 ✗ expansion_nt expansion_nt::operator+ (double rhs) const {
137 expansion* e = expansion::new_expansion_on_heap(
138 expansion::sum_capacity(rep(), rhs)
139 );
140 ✗ e->assign_sum(rep(), rhs);
141 ✗ return expansion_nt(e);
142 }
143
144 ✗ expansion_nt expansion_nt::operator- (double rhs) const {
145 expansion* e = expansion::new_expansion_on_heap(
146 expansion::diff_capacity(rep(), rhs)
147 );
148 ✗ e->assign_diff(rep(), rhs);
149 ✗ return expansion_nt(e);
150 }
151
152 ✗ expansion_nt expansion_nt::operator* (double rhs) const {
153 expansion* e = expansion::new_expansion_on_heap(
154 expansion::product_capacity(rep(), rhs)
155 );
156 ✗ e->assign_product(rep(), rhs);
157 ✗ return expansion_nt(e);
158 }
159
160 /************************************************************************/
161
162 67558 expansion_nt expansion_nt::operator- () const {
163 expansion_nt result(*this);
164 result.rep().negate();
165 67558 return result;
166 }
167
168 /************************************************************************/
169
170 2441004 expansion_nt expansion_nt_determinant(
171 const expansion_nt& a00,const expansion_nt& a01,
172 const expansion_nt& a10,const expansion_nt& a11
173 ) {
174 expansion* result = expansion::new_expansion_on_heap(
175 expansion::det2x2_capacity(a00.rep(),a01.rep(),a10.rep(),a11.rep())
176 );
177 2441004 result->assign_det2x2(a00.rep(),a01.rep(),a10.rep(),a11.rep());
178 2441004 return expansion_nt(result);
179 }
180
181
182 17274 expansion_nt expansion_nt_determinant(
183 const expansion_nt& a00,const expansion_nt& a01,const expansion_nt& a02,
184 const expansion_nt& a10,const expansion_nt& a11,const expansion_nt& a12,
185 const expansion_nt& a20,const expansion_nt& a21,const expansion_nt& a22
186 ) {
187 // First compute the det2x2
188 const expansion& m01 =
189 17274 expansion_det2x2(a00.rep(), a10.rep(), a01.rep(), a11.rep());
190 const expansion& m02 =
191 17274 expansion_det2x2(a00.rep(), a20.rep(), a01.rep(), a21.rep());
192 const expansion& m12 =
193 17274 expansion_det2x2(a10.rep(), a20.rep(), a11.rep(), a21.rep());
194
195 // Now compute the minors of rank 3
196 17274 const expansion& z1 = expansion_product(m01,a22.rep());
197 17274 const expansion& z2 = expansion_product(m02,a12.rep()).negate();
198 17274 const expansion& z3 = expansion_product(m12,a02.rep());
199
200 17274 return expansion_nt(expansion_nt::SUM, z1, z2, z3);
201 }
202
203 ✗ expansion_nt expansion_nt_determinant(
204 const expansion_nt& a00,const expansion_nt& a01,
205 const expansion_nt& a02,const expansion_nt& a03,
206 const expansion_nt& a10,const expansion_nt& a11,
207 const expansion_nt& a12,const expansion_nt& a13,
208 const expansion_nt& a20,const expansion_nt& a21,
209 const expansion_nt& a22,const expansion_nt& a23,
210 const expansion_nt& a30,const expansion_nt& a31,
211 const expansion_nt& a32,const expansion_nt& a33
212 ) {
213
214 // First compute the det2x2
215 const expansion& m01 =
216 ✗ expansion_det2x2(a10.rep(),a00.rep(),a11.rep(),a01.rep());
217 const expansion& m02 =
218 ✗ expansion_det2x2(a20.rep(),a00.rep(),a21.rep(),a01.rep());
219 const expansion& m03 =
220 ✗ expansion_det2x2(a30.rep(),a00.rep(),a31.rep(),a01.rep());
221 const expansion& m12 =
222 ✗ expansion_det2x2(a20.rep(),a10.rep(),a21.rep(),a11.rep());
223 const expansion& m13 =
224 ✗ expansion_det2x2(a30.rep(),a10.rep(),a31.rep(),a11.rep());
225 const expansion& m23 =
226 ✗ expansion_det2x2(a30.rep(),a20.rep(),a31.rep(),a21.rep());
227
228 // Now compute the minors of rank 3
229 ✗ const expansion& m012_1 = expansion_product(m12,a02.rep());
230 ✗ expansion& m012_2 = expansion_product(m02,a12.rep()); m012_2.negate();
231 ✗ const expansion& m012_3 = expansion_product(m01,a22.rep());
232 ✗ const expansion& m012 = expansion_sum3(m012_1, m012_2, m012_3);
233
234 ✗ const expansion& m013_1 = expansion_product(m13,a02.rep());
235 ✗ expansion& m013_2 = expansion_product(m03,a12.rep()); m013_2.negate();
236
237 ✗ const expansion& m013_3 = expansion_product(m01,a32.rep());
238 ✗ const expansion& m013 = expansion_sum3(m013_1, m013_2, m013_3);
239
240 ✗ const expansion& m023_1 = expansion_product(m23,a02.rep());
241 ✗ expansion& m023_2 = expansion_product(m03,a22.rep()); m023_2.negate();
242 ✗ const expansion& m023_3 = expansion_product(m02,a32.rep());
243 ✗ const expansion& m023 = expansion_sum3(m023_1, m023_2, m023_3);
244
245 ✗ const expansion& m123_1 = expansion_product(m23,a12.rep());
246 ✗ expansion& m123_2 = expansion_product(m13,a22.rep()); m123_2.negate();
247 ✗ const expansion& m123_3 = expansion_product(m12,a32.rep());
248 ✗ const expansion& m123 = expansion_sum3(m123_1, m123_2, m123_3);
249
250 // Now compute the minors of rank 4
251 ✗ const expansion& m0123_1 = expansion_product(m123,a03.rep());
252 ✗ const expansion& m0123_2 = expansion_product(m023,a13.rep());
253 ✗ const expansion& m0123_3 = expansion_product(m013,a23.rep());
254 ✗ const expansion& m0123_4 = expansion_product(m012,a33.rep());
255
256 ✗ const expansion& z1 = expansion_sum(m0123_1, m0123_3);
257 ✗ const expansion& z2 = expansion_sum(m0123_2, m0123_4);
258
259 ✗ return expansion_nt(expansion_nt::DIFF,z1,z2);
260 }
261
262 /***********************************************************************/
263
264 namespace Numeric {
265
266
1/2
✓ Branch 0 taken 2235047 times.
✗ Branch 1 not taken.
 2235047 template<> Sign ratio_compare(
267 const expansion_nt& a_num, const expansion_nt& a_denom,
268 const expansion_nt& b_num, const expansion_nt& b_denom
269 ) {
270 // TODO HERE: CHECK THAT THIS FITS ON STACK
271
272
1/2
✓ Branch 0 taken 2235047 times.
✗ Branch 1 not taken.
 2235047 Sign s1 = Sign(a_num.sign()*a_denom.sign());
273 2235047 Sign s2 = Sign(b_num.sign()*b_denom.sign());
274
2/2
✓ Branch 0 taken 2219493 times.
✓ Branch 1 taken 15554 times.
 2235047 if(s1 == ZERO && s2 == ZERO) {
275 return ZERO;
276 }
277
2/2
✓ Branch 0 taken 276414 times.
✓ Branch 1 taken 1943079 times.
 2219493 if(s1 != s2) {
278
2/2
✓ Branch 0 taken 162487 times.
✓ Branch 1 taken 113927 times.
 438901 return (int(s1) > int(s2) ? POSITIVE : NEGATIVE);
279 }
280
2/2
✓ Branch 0 taken 239043 times.
✓ Branch 1 taken 1704036 times.
 1943079 if(a_denom == b_denom) {
281 239043 const expansion& diff_num = expansion_diff(
282 a_num.rep(), b_num.rep()
283 );
284 239043 return Sign(diff_num.sign() * a_denom.sign());
285 }
286 1704036 const expansion& num_a = expansion_product(
287 a_num.rep(), b_denom.rep()
288 );
289 1704036 const expansion& num_b = expansion_product(
290 b_num.rep(), a_denom.rep()
291 );
292 1704036 const expansion& diff_num = expansion_diff(num_a, num_b);
293 return Sign(
294
1/2
✓ Branch 0 taken 1704036 times.
✗ Branch 1 not taken.
 1704036 diff_num.sign() * a_denom.sign() * b_denom.sign()
295 1704036 );
296 }
297
298 }
299
300 }
301