Geogram  Version 1.9.1-rc
A programming library of geometric algorithms
expansion_nt.h
Go to the documentation of this file.
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 #ifndef GEOGRAM_NUMERICS_EXPANSION_NT
41 #define GEOGRAM_NUMERICS_EXPANSION_NT
42 
43 #include <geogram/basic/common.h>
45 #include <geogram/basic/matrix.h>
47 
56 namespace GEO {
57 
58  class expansion_nt;
59 
70  class GEOGRAM_API expansion_nt {
71  public:
72 
77  enum Operation {
78  SUM, DIFF, PRODUCT
79  };
80 
84  expansion_nt() : rep_(nullptr) {
85  }
86 
91  explicit expansion_nt(double x) {
93  rep()[0] = x;
94  rep().set_length(1);
95  }
96 
103  explicit expansion_nt(const expansion& rhs) {
105  rep().assign(rhs);
106  }
107 
121  explicit expansion_nt(
122  Operation op, const expansion& x, const expansion& y
123  ) {
124  switch(op) {
125  case SUM:
128  );
129  rep_->assign_sum(x,y);
130  break;
131  case DIFF:
134  );
135  rep_->assign_diff(x,y);
136  break;
137  case PRODUCT:
140  );
141  rep_->assign_product(x,y);
142  break;
143  }
144  }
145 
159  explicit expansion_nt(
160  Operation op,
161  const expansion& x, const expansion& y, const expansion& z
162  ) {
163  switch(op) {
164  case SUM:
167  );
168  rep_->assign_sum(x,y,z);
169  break;
170  case DIFF:
172  break;
173  case PRODUCT:
176  );
177  rep_->assign_product(x,y,z);
178  break;
179  }
180  }
181 
195  explicit expansion_nt(
196  Operation op,
197  const expansion& x, const expansion& y,
198  const expansion& z, const expansion& t
199  ) {
200  switch(op) {
201  case SUM:
203  expansion::sum_capacity(x,y,z,t)
204  );
205  rep_->assign_sum(x,y,z,t);
206  break;
207  case DIFF:
209  break;
210  case PRODUCT:
211  // HERE: TODO CHECK SIZE
212  const expansion& p1 = expansion_product(x,y);
213  const expansion& p2 = expansion_product(z,t);
216  );
217  rep_->assign_sum(p1,p2);
218  break;
219  }
220  }
221 
235  explicit expansion_nt(Operation op, double x, double y) {
236  switch(op) {
237  case SUM:
240  );
241  rep_->assign_sum(x,y);
242  break;
243  case DIFF:
246  );
247  rep_->assign_diff(x,y);
248  break;
249  case PRODUCT:
252  );
253  rep_->assign_product(x,y);
254  break;
255  }
256  }
257 
263  copy(rhs);
264  }
265 
272  rep_ = nullptr;
273  std::swap(rep_, rhs.rep_);
274  }
275 
281  expansion_nt& operator= (const expansion_nt& rhs) {
282  if(&rhs != this) {
283  cleanup();
284  copy(rhs);
285  }
286  return *this;
287  }
288 
294  expansion_nt& operator= (expansion_nt&& rhs) {
295  if(&rhs != this) {
296  cleanup();
297  std::swap(rep_, rhs.rep_);
298  }
299  return *this;
300  }
301 
308  cleanup();
309  }
310 
316  void optimize() {
317  rep().optimize();
318  }
319 
323  void negate() {
324  rep().negate();
325  }
326 
327  /********************************************************************/
328 
334  expansion_nt& operator+= (const expansion_nt& rhs);
335 
341  expansion_nt& operator-= (const expansion_nt& rhs);
342 
348  expansion_nt& operator*= (const expansion_nt& rhs);
349 
355  expansion_nt& operator+= (double rhs);
356 
362  expansion_nt& operator-= (double rhs);
363 
372  expansion_nt& operator*= (double rhs);
373 
374  /********************************************************************/
375 
382 
390 
397  expansion_nt operator* (const expansion_nt& rhs) const;
398 
404  expansion_nt operator+ (double rhs) const;
405 
411  expansion_nt operator- (double rhs) const;
412 
418  expansion_nt operator* (double rhs) const;
419 
420  /********************************************************************/
421 
427 
428  /********************************************************************/
429 
434  Sign compare(const expansion_nt& rhs) const {
435  return rep().compare(rhs.rep());
436  }
437 
442  Sign compare(double rhs) const {
443  return rep().compare(rhs);
444  }
445 
453  bool operator> (const expansion_nt& rhs) const {
454  return (int(compare(rhs))>0);
455  }
456 
464  bool operator>= (const expansion_nt& rhs) const {
465  return (int(compare(rhs))>=0);
466  }
467 
475  bool operator< (const expansion_nt& rhs) const {
476  return (int(compare(rhs))<0);
477  }
478 
486  bool operator<= (const expansion_nt& rhs) const {
487  return (int(compare(rhs))<=0);
488  }
489 
497  bool operator> (double rhs) const {
498  return (int(compare(rhs))>0);
499  }
500 
508  bool operator>= (double rhs) const {
509  return (int(compare(rhs))>=0);
510  }
511 
519  bool operator< (double rhs) const {
520  return (int(compare(rhs))<0);
521  }
522 
530  bool operator<= (double rhs) const {
531  return (int(compare(rhs))<=0);
532  }
533 
534  /********************************************************************/
535 
541  double estimate() const {
542  return rep().estimate();
543  }
544 
549  Sign sign() const {
550  return rep().sign();
551  }
552 
560  index_t length() const {
561  return rep().length();
562  }
563 
572  double component(index_t i) const {
573  geo_debug_assert(i < length());
574  return rep()[i];
575  }
576 
586  rep_(rep) {
587  }
588 
598  return *rep_;
599  }
600 
609  const expansion& rep() const {
610  return *rep_;
611  }
612 
618  std::string to_string() const {
619  return (rep_ == nullptr) ?
620  std::string("null") :
621  rep_->to_string() ;
622  }
623 
624  protected:
625 
632  void copy(const expansion_nt& rhs) {
633  if(rhs.rep_ == nullptr) {
634  rep_ = nullptr;
635  } else {
637  rep_->set_length(rhs.rep().length());
638  for(index_t i=0; i<rep_->length(); ++i) {
639  (*rep_)[i] = rhs.rep()[i];
640  }
641  }
642  }
643 
647  void cleanup() {
648  if(rep_ != nullptr) {
650  rep_ = nullptr;
651  }
652  }
653 
654  private:
655  expansion* rep_;
656  friend expansion_nt operator- (double a, const expansion_nt& b);
657 
658  friend expansion_nt expansion_nt_sq_dist(
659  const double* a, const double* b, coord_index_t dim
660  );
661 
662  friend expansion_nt expansion_nt_dot_at(
663  const double* a, const double* b, const double* c,
664  coord_index_t dim
665  );
666  // friend class rational_nt;
667  };
668 
676  inline expansion_nt operator+ (double a, const expansion_nt& b) {
677  return b + a;
678  }
679 
687  inline expansion_nt operator- (double a, const expansion_nt& b) {
688  expansion_nt result = b - a;
689  result.rep().negate();
690  return result;
691  }
692 
700  inline expansion_nt operator* (double a, const expansion_nt& b) {
701  return b * a;
702  }
703 
712  inline bool operator== (const expansion_nt& a, const expansion_nt& b) {
713  return a.rep().equals(b.rep());
714  }
715 
724  inline bool operator== (const expansion_nt& a, double b) {
725  return a.rep().equals(b);
726  }
727 
736  inline bool operator== (double a, const expansion_nt& b) {
737  return b.rep().equals(a);
738  }
739 
748  inline bool operator!= (const expansion_nt& a, const expansion_nt& b) {
749  return !a.rep().equals(b.rep());
750  }
751 
760  inline bool operator!= (const expansion_nt& a, double b) {
761  return !a.rep().equals(b);
762  }
763 
772  inline bool operator!= (double a, const expansion_nt& b) {
773  return !b.rep().equals(a);
774  }
775 
787  const double* a, const double* b, coord_index_t dim
788  ) {
791  );
792  result->assign_sq_dist(a, b, dim);
793  return expansion_nt(result);
794  }
795 
808  const double* a, const double* b, const double* c, coord_index_t dim
809  ) {
812  );
813  result->assign_dot_at(a, b, c, dim);
814  return expansion_nt(result);
815  }
816 
817  /************************************************************************/
818 
824  template <> inline Sign geo_sgn(const expansion_nt& x) {
825  return x.sign();
826  }
827 
835  template <> inline Sign geo_cmp(
836  const expansion_nt& x, const expansion_nt& y
837  ) {
838  return x.compare(y);
839  }
840 
841  /************************************************************************/
842 
850  inline bool expansion_nt_is_zero(const expansion_nt& x) {
851  return (x.sign() == GEO::ZERO);
852  }
853 
861  inline bool expansion_nt_is_one(const expansion_nt& x) {
862  return x.rep().equals(1.0);
863  }
864 
865 
875  const expansion_nt& x, const expansion_nt& y
876  ) {
877  const expansion& diff = expansion_diff(x.rep(), y.rep());
878  return diff.sign();
879  }
880 
888  expansion_nt result(
891  ))
892  );
893  result.rep().assign_square(x.rep());
894  return result;
895  }
896 
897 
905  const expansion_nt& a00,const expansion_nt& a01,
906  const expansion_nt& a10,const expansion_nt& a11
907  );
908 
916  const expansion_nt& a00,const expansion_nt& a01,const expansion_nt& a02,
917  const expansion_nt& a10,const expansion_nt& a11,const expansion_nt& a12,
918  const expansion_nt& a20,const expansion_nt& a21,const expansion_nt& a22
919  );
920 
928  const expansion_nt& a00,const expansion_nt& a01,
929  const expansion_nt& a02,const expansion_nt& a03,
930  const expansion_nt& a10,const expansion_nt& a11,
931  const expansion_nt& a12,const expansion_nt& a13,
932  const expansion_nt& a20,const expansion_nt& a21,
933  const expansion_nt& a22,const expansion_nt& a23,
934  const expansion_nt& a30,const expansion_nt& a31,
935  const expansion_nt& a32,const expansion_nt& a33
936  );
937 
938 // Make things a bit faster if target OS has large stack size
939 #ifdef GEO_HAS_BIG_STACK
940 
946  template <> inline expansion_nt det2x2(
947  const expansion_nt& a11, const expansion_nt& a12,
948  const expansion_nt& a21, const expansion_nt& a22
949  ) {
951  a11,a12,
952  a21,a22
953  );
954  }
955 
961  template <> inline expansion_nt det3x3(
962  const expansion_nt& a11, const expansion_nt& a12,
963  const expansion_nt& a13,
964  const expansion_nt& a21, const expansion_nt& a22,
965  const expansion_nt& a23,
966  const expansion_nt& a31, const expansion_nt& a32,
967  const expansion_nt& a33
968  ) {
970  a11,a12,a13,
971  a21,a22,a23,
972  a31,a32,a33
973  );
974  }
975 
981  template <> inline expansion_nt det4x4(
982  const expansion_nt& a11, const expansion_nt& a12,
983  const expansion_nt& a13, const expansion_nt& a14,
984  const expansion_nt& a21, const expansion_nt& a22,
985  const expansion_nt& a23, const expansion_nt& a24,
986  const expansion_nt& a31, const expansion_nt& a32,
987  const expansion_nt& a33, const expansion_nt& a34,
988  const expansion_nt& a41, const expansion_nt& a42,
989  const expansion_nt& a43, const expansion_nt& a44
990  ) {
992  a11,a12,a13,a14,
993  a21,a22,a23,a24,
994  a31,a32,a33,a34,
995  a41,a42,a43,a44
996  );
997  }
998 
999 #endif
1000 
1001  /************************************************************************/
1002 }
1003 
1010 inline std::ostream& operator<< (
1011  std::ostream& os, const GEO::expansion_nt& a
1012 ) {
1013  return os << a.estimate();
1014 }
1015 
1023 inline std::istream& operator>> ( std::istream& is, GEO::expansion_nt& a) {
1024  double d;
1025  is >> d;
1026  if (is) {
1027  a = GEO::expansion_nt(d);
1028  }
1029  return is;
1030 }
1031 
1032 /*****************************************************************************/
1033 
1034 namespace GEO {
1035 
1036  /**************************************************************************/
1037 
1038  namespace Numeric {
1039 
1040  template<> inline void optimize_number_representation(expansion_nt& x) {
1041  x.optimize();
1042  }
1043 
1051  template<> Sign GEOGRAM_API ratio_compare(
1052  const expansion_nt& a_num, const expansion_nt& a_denom,
1053  const expansion_nt& b_num, const expansion_nt& b_denom
1054  );
1055  }
1056 
1057  /**************************************************************************/
1058 
1060 }
1061 
1062 #endif
#define geo_assert_not_reached
Sets a non reachable point in the program.
Definition: assert.h:177
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition: assert.h:196
std::ostream & operator<<(std::ostream &output, const Matrix< DIM, FT > &m)
Writes a matrix to a stream.
Definition: matrix.h:466
std::istream & operator>>(std::istream &input, Matrix< DIM, FT > &m)
Reads a matrix from a stream.
Definition: matrix.h:489
Expansion_nt (expansion Number Type) is used to compute the sign of polynoms exactly.
Definition: expansion_nt.h:70
const expansion & rep() const
Gets the internal expansion that represents this expansion_nt.
Definition: expansion_nt.h:609
double estimate() const
Computes an approximation of the stored value in this expansion.
Definition: expansion_nt.h:541
void copy(const expansion_nt &rhs)
Copies an expansion into this one.
Definition: expansion_nt.h:632
Sign sign() const
Gets the sign of this expansion_nt.
Definition: expansion_nt.h:549
expansion_nt expansion_nt_sq_dist(const double *a, const double *b, coord_index_t dim)
Computes an expansion that represents the square distance between two points.
Definition: expansion_nt.h:786
double component(index_t i) const
Gets the i-th component of this expansion.
Definition: expansion_nt.h:572
~expansion_nt()
Expansion_nt destructor.
Definition: expansion_nt.h:307
std::string to_string() const
Gets a string representation of this expansion.
Definition: expansion_nt.h:618
expansion_nt expansion_nt_dot_at(const double *a, const double *b, const double *c, coord_index_t dim)
Computes an expansion that represents the dot product of two vectors determined by three points.
Definition: expansion_nt.h:807
expansion_nt(expansion_nt &&rhs)
Move-constructor.
Definition: expansion_nt.h:271
expansion_nt(Operation op, const expansion &x, const expansion &y, const expansion &z)
Constructs a new expansion_nt from three expansions.
Definition: expansion_nt.h:159
expansion_nt()
Constructs an uninitialized expansion_nt.
Definition: expansion_nt.h:84
expansion_nt(const expansion &rhs)
Constructs a new expansion_nt from an expansion.
Definition: expansion_nt.h:103
void cleanup()
Cleanups the memory associated with this expansion_nt.
Definition: expansion_nt.h:647
index_t length() const
Gets the length of this expansion.
Definition: expansion_nt.h:560
void optimize()
Optimizes the internal representation without changing the represented value.
Definition: expansion_nt.h:316
expansion_nt(Operation op, const expansion &x, const expansion &y)
Constructs a new expansion_nt from two expansions.
Definition: expansion_nt.h:121
Sign compare(const expansion_nt &rhs) const
Compares two expansion_nt.
Definition: expansion_nt.h:434
expansion_nt(expansion *rep)
Constructs a new expansion_nt from an expansion.
Definition: expansion_nt.h:585
expansion_nt(Operation op, double x, double y)
Constructs a new expansion_nt from two doubles.
Definition: expansion_nt.h:235
expansion_nt(Operation op, const expansion &x, const expansion &y, const expansion &z, const expansion &t)
Constructs a new expansion_nt from four expansions.
Definition: expansion_nt.h:195
Operation
This type is used by the constructor that takes two expansions.
Definition: expansion_nt.h:77
expansion_nt(const expansion_nt &rhs)
Copy-constructor.
Definition: expansion_nt.h:262
Sign compare(double rhs) const
Compares an expansion_nt with a double.
Definition: expansion_nt.h:442
void negate()
Flips the sign of an expansion.
Definition: expansion_nt.h:323
expansion & rep()
Gets the internal expansion that represents this expansion_nt.
Definition: expansion_nt.h:597
expansion_nt(double x)
Constructs a new expansion_nt from a double.
Definition: expansion_nt.h:91
Represents numbers in arbitrary precision with a low-level API.
expansion & assign_sq_dist(const double *p1, const double *p2, coord_index_t dim)
Assigns the squared distance between two points to this expansion (should not be used by client code)...
static index_t product_capacity(double a, double b)
Computes the required capacity of an expansion to store the exact product of two doubles.
expansion & negate()
Changes the sign of an expansion.
expansion & assign_dot_at(const double *p1, const double *p2, const double *p0, coord_index_t dim)
Assigns the dot product of two vectors to this expansion (should not be used by client code).
static index_t sum_capacity(double a, double b)
Computes the required capacity to store the sum of two doubles.
index_t length() const
Gets the length of this expansion.
static void delete_expansion_on_heap(expansion *e)
Deallocates an expansion on the heap.
static index_t diff_capacity(double a, double b)
Computes the required capacity of an expansion to store the exact difference of two doubles.
bool equals(const expansion &rhs) const
Compares two expansions.
index_t capacity() const
Gets the capacity of this expansion.
static index_t sq_dist_capacity(coord_index_t dim)
Computes the required capacity of an expansion to store the exact squared distance between two points...
static expansion * new_expansion_on_heap(index_t capa)
Allocates an expansion on the heap.
static index_t square_capacity(double a)
Computes the required capacity of an expansion to store the exact square of a double.
static index_t dot_at_capacity(coord_index_t dim)
Computes the required capacity of an expansion to store the exact dot product between two vectors.
expansion & assign_square(double a)
Assigns the square of a double to this expansion (should not be used by client code).
Sign sign() const
Gets the sign of the expansion.
rationalg (generic rational) is used to compute the sign of rational fractions exactly.
Definition: rationalg.h:59
Common include file, providing basic definitions. Should be included before anything else by all head...
Generic matrix type.
Implementation of multi-precision arithmetics.
bool operator!=(const aligned_allocator< T1, A1 > &, const aligned_allocator< T2, A2 > &)
Tests whether two aligned_allocators are different.
Definition: memory.h:617
bool operator==(const aligned_allocator< T1, A1 > &, const aligned_allocator< T2, A2 > &)
Tests whether two aligned_allocators are equal.
Definition: memory.h:606
void copy(void *to, const void *from, size_t size)
Copies a memory block.
Definition: memory.h:129
Sign ratio_compare(const T &a_num, const T &a_denom, const T &b_num, const T &b_denom)
Compares two rational numbers given as separate numerators and denominators.
Definition: numeric.h:278
void optimize_number_representation(T &x)
place holder for optimizing internal number representation
Definition: numeric.h:267
Global Vorpaline namespace.
Definition: basic.h:55
bool expansion_nt_is_zero(const expansion_nt &x)
Tests whether an expansion_nt is zero.
Definition: expansion_nt.h:850
Quaternion operator-(const Quaternion &a, const Quaternion &b)
Computes the difference between two Quaternion.
Definition: quaternion.h:252
expansion_nt expansion_nt_determinant(const expansion_nt &a00, const expansion_nt &a01, const expansion_nt &a10, const expansion_nt &a11)
Computes a 2x2 determinant.
T det3x3(const T &a11, const T &a12, const T &a13, const T &a21, const T &a22, const T &a23, const T &a31, const T &a32, const T &a33)
Computes a three-by-three determinant.
Definition: determinant.h:69
T det4x4(const T &a11, const T &a12, const T &a13, const T &a14, const T &a21, const T &a22, const T &a23, const T &a24, const T &a31, const T &a32, const T &a33, const T &a34, const T &a41, const T &a42, const T &a43, const T &a44)
Computes a four-by-four determinant.
Definition: determinant.h:85
Sign geo_sgn(const T &x)
Gets the sign of a value.
Definition: numeric.h:90
expansion_nt expansion_nt_square(const expansion_nt &x)
Computes the square of an expansion_nt.
Definition: expansion_nt.h:887
geo_index_t index_t
The type for storing and manipulating indices.
Definition: numeric.h:329
Sign expansion_nt_compare(const expansion_nt &x, const expansion_nt &y)
Compares two expansion_nt.
Definition: expansion_nt.h:874
Sign
Integer constants that represent the sign of a value.
Definition: numeric.h:68
@ ZERO
Definition: numeric.h:72
T det2x2(const T &a11, const T &a12, const T &a21, const T &a22)
Computes a two-by-two determinant.
Definition: determinant.h:58
Sign geo_cmp(const T &a, const T &b)
Compares two values.
Definition: numeric.h:106
bool expansion_nt_is_one(const expansion_nt &x)
Tests whether an expansion_nt is equal to one.
Definition: expansion_nt.h:861
Quaternion operator+(const Quaternion &a, const Quaternion &b)
Computes the sum of two Quaternion.
Definition: quaternion.h:239
geo_coord_index_t coord_index_t
The type for storing coordinate indices, and iterating on the coordinates of a point.
Definition: numeric.h:363
Generic implementation of rational type.