Geogram  Version 1.9.1-rc
A programming library of geometric algorithms
interval_nt.h
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 GEO_INTERVAL_NT
41 #define GEO_INTERVAL_NT
42 
44 #include <geogram/basic/vechg.h>
45 #include <iomanip>
46 #include <limits>
47 #include <cmath>
48 
49 #ifdef __SSE__
50 #include <xmmintrin.h>
51 #else
52 #include <fenv.h>
53 #endif
54 
55 // Uncomment to activate checks (keeps an arbitrary precision
56 // representation of the number and checks that the interval
57 // contains it).
58 
59 // #define INTERVAL_CHECK
60 
61 namespace GEO {
62 
71  class intervalBase {
72  public:
73 
74  static void set_FPU_round_to_nearest() {
75 #ifdef __SSE__
76  _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST);
77 #else
78  fesetround(FE_TONEAREST);
79 #endif
80  }
81 
82  static void set_FPU_round_to_upper() {
83 #ifdef __SSE__
84  _MM_SET_ROUNDING_MODE(_MM_ROUND_UP);
85 #else
86  fesetround(FE_UPWARD);
87 #endif
88  }
89 
90  enum Sign2 {
91  SIGN2_ERROR = -1,
92  SIGN2_ZERO = 0,
93  SIGN2_NP,
94  SIGN2_PP,
95  SIGN2_ZP,
96  SIGN2_NN,
97  SIGN2_NZ
98  };
99 
100  static bool sign_is_determined(Sign2 s) {
101  return
102  s == SIGN2_ZERO ||
103  s == SIGN2_NN ||
104  s == SIGN2_PP ;
105  }
106 
107  static bool sign_is_non_zero(Sign2 s) {
108  return
109  s == SIGN2_NN ||
110  s == SIGN2_PP ;
111  }
112 
113  static Sign convert_sign(Sign2 s) {
114  geo_assert(sign_is_determined(s));
115  if(s == SIGN2_NN) {
116  return NEGATIVE;
117  }
118  if(s == SIGN2_PP) {
119  return POSITIVE;
120  }
121  return ZERO;
122  }
123 
124  intervalBase() {
125  control_set(0);
126  }
127 
128  intervalBase(double x) {
129  control_set(x);
130  }
131 
132  intervalBase(const intervalBase& rhs) = default;
133 
134  intervalBase& operator=(const intervalBase& rhs) = default;
135 
136  protected:
137 #ifdef INTERVAL_CHECK
138  void control_set(const expansion_nt& x) {
139  control_ = x;
140  }
141  void control_set(const intervalBase& x) {
142  control_ = x.control_;
143  }
144  void control_set(double x) {
145  control_ = expansion_nt(x);
146  }
147  void control_negate() {
148  control_.rep().negate();
149  }
150  void control_add(const intervalBase& x) {
151  control_ += x.control_;
152  }
153  void control_sub(const intervalBase& x) {
154  control_ -= x.control_;
155  }
156  void control_mul(const intervalBase& x) {
157  control_ *= x.control_;
158  }
159  void control_check(double inf, double sup) {
160  typedef std::numeric_limits< double > dbl;
161  if(inf > sup) {
162  std::cerr.precision(dbl::max_digits10);
163  std::cerr << "inf() > sup() !!" << std::endl;
164  std::cerr << "inf()=" << inf << std::endl;
165  std::cerr << "sup()=" << sup << std::endl;
167  }
168  if(control_ < inf || control_ > sup) {
169  std::cerr.precision(dbl::max_digits10);
170  std::cerr << "[" << inf << "," << sup << "]"
171  << " " << control_.estimate() << ":"
172  << control_.rep().length()
173  << std::endl;
175  }
176  }
177  expansion_nt control_;
178 #else
179  void control_set(double x) {
180  geo_argused(x);
181  }
182  void control_set(const expansion_nt& x) {
183  geo_argused(x);
184  }
185  void control_set(const intervalBase& x) {
186  geo_argused(x);
187  }
188  void control_negate() {
189  }
190  void control_add(const intervalBase& x) {
191  geo_argused(x);
192  }
193  void control_sub(const intervalBase& x) {
194  geo_argused(x);
195  }
196  void control_mul(const intervalBase& x) {
197  geo_argused(x);
198  }
199  void control_check(double inf, double sup) {
200  geo_argused(inf);
201  geo_argused(sup);
202  }
203 #endif
204  };
205 
206 
207 /*******************************************************************/
208 
212  class intervalRU : public intervalBase {
213  public:
214 
227  struct Rounding {
228  Rounding() {
229  set_FPU_round_to_upper();
230  }
231  ~Rounding() {
232  set_FPU_round_to_nearest();
233  }
234  };
235 
236  intervalRU() :
237  intervalBase(),
238  ln_(0.0),
239  u_(0.0)
240  {
241  control_check();
242  }
243 
244  intervalRU(double x) :
245  intervalBase(x),
246  ln_(-x),
247  u_(x)
248  {
249  control_check();
250  }
251 
252  intervalRU(double l, double u) : ln_(-l), u_(u) {
253  // note: we cannot control check here.
254  }
255 
256  intervalRU(const intervalRU& rhs) = default;
257 
258  intervalRU(const expansion_nt& rhs) {
259  *this = rhs;
260  }
261 
262  intervalRU& operator=(const intervalRU& rhs) = default;
263 
264  intervalRU& operator=(double rhs) {
265  ln_ = -rhs;
266  u_ = rhs;
267  control_set(rhs);
268  control_check();
269  return *this;
270  }
271 
272  intervalRU& operator=(const expansion_nt& rhs) {
273 
274  // Optimized expansion-to-interval conversion:
275  //
276  // Add components starting from the one of largest magnitude
277  // Stop as soon as next component is smaller than ulp (and then
278  // expand interval by ulp).
279 
280  index_t l = rhs.length();
281  ln_ = -rhs.component(l-1);
282  u_ = rhs.component(l-1);
283 
284  for(int comp_idx=int(l)-2; comp_idx>=0; --comp_idx) {
285  double comp = rhs.component(index_t(comp_idx));
286  u_ += comp;
287  ln_ -= comp;
288  if(comp > 0) {
289  double new_u = u_ + comp;
290  if(new_u == u_) {
291  u_ = std::nextafter(
292  u_, std::numeric_limits<double>::infinity()
293  );
294  break;
295  } else {
296  u_ = new_u;
297  }
298  } else {
299  // If we stored l, we would write:
300  // new_l = l + comp
301  // But we store ln = -l, so we write:
302  // -new_ln = -ln + comp
303  // Which means:
304  // new_ln = ln - comp
305 
306  double new_ln = ln_ - comp;
307  if(new_ln == ln_) {
308  ln_ = std::nextafter(
309  ln_, std::numeric_limits<double>::infinity()
310  );
311  break;
312  } else {
313  ln_ = new_ln;
314  }
315  }
316  }
317 
318  control_set(rhs);
319  control_check();
320  return *this;
321  }
322 
323  double inf() const {
324  return -ln_;
325  }
326 
327  double sup() const {
328  return u_;
329  }
330 
331  double estimate() const {
332  // 0.5*(lb+ub) ->
333  return 0.5*(-ln_+u_);
334  }
335 
336  bool is_nan() const {
337  return !(ln_==ln_) || !(u_==u_);
338  }
339 
340  Sign2 sign() const {
341  // Branchless
342  int lz = int(ln_ == 0);
343  int ln = int(ln_ > 0); // inverted, it is ln_ !!!
344  int lp = int(ln_ < 0); // inverted, it is ln_ !!!
345  int uz = int(u_ == 0);
346  int un = int(u_ < 0);
347  int up = int(u_ > 0);
348  Sign2 result = Sign2(
349  ln*up*SIGN2_NP+
350  lp*up*SIGN2_PP+
351  lz*up*SIGN2_ZP+
352  ln*un*SIGN2_NN+
353  ln*uz*SIGN2_NZ
354  );
355  result = Sign2(
356  int(result) +
357  int(result==SIGN2_ZERO && !(lz&&uz)) * SIGN2_ERROR
358  );
359  return result;
360  }
361 
362  intervalRU& negate() {
363  std::swap(ln_, u_);
364  control_negate();
365  control_check();
366  return *this;
367  }
368 
369  intervalRU& operator+=(const intervalRU &x) {
370  // lb += x.lb -> -lbn += -x.lbn -> lbn += x.lbn
371  ln_ += x.ln_;
372  u_ += x.u_;
373  control_add(x);
374  control_check();
375  return *this;
376  }
377 
378  intervalRU& operator-=(const intervalRU &x) {
379  // +=(x.negate()) ->
380  ln_ += x.u_;
381  u_ += x.ln_;
382  control_sub(x);
383  control_check();
384  return *this;
385  }
386 
387 
388  intervalRU& operator*=(const intervalRU &b) {
389  // get bounds of both operands
390  double aln = ln_;
391  double au = u_;
392  double bln = b.ln_;
393  double bu = b.u_;
394 
395  // negated bounds round to upper
396  // (equivalent to bounds round to lower)
397  double lln = (-aln)*bln;
398  double lun = aln*bu;
399  double uln = au*bln;
400  double uun = (-au)*bu;
401 
402  // bounds round to upper
403  double ll = aln*bln;
404  double lu = (-aln)*bu;
405  double ul = au*(-bln);
406  double uu = au*bu;
407 
408  ln_ = std::max(std::max(lln,lun),std::max(uln,uun));
409  u_ = std::max(std::max(ll,lu),std::max(ul,uu));
410 
411  control_mul(b);
412  control_check();
413  return *this;
414  }
415 
416  protected:
417 
418 #ifdef INTERVAL_CHECK
419  void control_check() {
420  // expansion_nt used in control_check() operates
421  // in round to nearest mode !!
422  set_FPU_round_to_nearest();
423  intervalBase::control_check(inf(),sup());
424  set_FPU_round_to_upper();
425  }
426 #else
427  void control_check() {
428  }
429 #endif
430  double ln_;
431  double u_;
432  };
433 
434 
435  inline intervalRU operator+(const intervalRU& a, const intervalRU& b) {
436  intervalRU result = a;
437  return result += b;
438  }
439 
440  inline intervalRU operator-(const intervalRU& a, const intervalRU& b) {
441  intervalRU result = a;
442  return result -= b;
443  }
444 
445  inline intervalRU operator*(const intervalRU& a, const intervalRU& b) {
446  intervalRU result = a;
447  return result *= b;
448  }
449 
450  /*************************************************************************/
451 
459  class intervalRN : public intervalBase {
460  public:
461 
462  // operates in default rounding mode
463  // (so Rounding subclass does nothing)
464  struct Rounding {
465  Rounding() {
466  }
467  ~Rounding() {
468  }
469  };
470 
471  intervalRN() :
472  intervalBase(),
473  lb_(0.0),
474  ub_(0.0)
475  {
476  control_check();
477  }
478 
479  intervalRN(double x) :
480  intervalBase(x),
481  lb_(x),
482  ub_(x)
483  {
484  control_check();
485  }
486 
487  intervalRN(double l, double u) : lb_(l), ub_(u) {
488  // note: we cannot control check here.
489  }
490 
491 
492  intervalRN(const intervalRN& rhs) = default;
493 
494  intervalRN(const expansion_nt& rhs) {
495  *this = rhs;
496  }
497 
498  intervalRN& operator=(const intervalRN& rhs) = default;
499 
500  intervalRN& operator=(double rhs) {
501  lb_ = rhs;
502  ub_ = rhs;
503  control_set(rhs);
504  control_check();
505  return *this;
506  }
507 
508  intervalRN& operator=(const expansion_nt& rhs) {
509 
510  // Optimized expansion-to-interval conversion:
511  //
512  // Add components starting from the one of largest magnitude
513  // Stop as soon as next component is smaller than ulp (and then
514  // expand interval by ulp).
515 
516  index_t l = rhs.length();
517  lb_ = rhs.component(l-1);
518  ub_ = rhs.component(l-1);
519 
520  for(int comp_idx=int(l)-2; comp_idx>=0; --comp_idx) {
521  double comp = rhs.component(index_t(comp_idx));
522  if(comp > 0) {
523  double nub = ub_ + comp;
524  if(nub == ub_) {
525  ub_ = std::nextafter(
526  ub_, std::numeric_limits<double>::infinity()
527  );
528  break;
529  } else {
530  ub_ = nub;
531  adjust();
532  }
533  } else {
534  double nlb = lb_ + comp;
535  if(nlb == lb_) {
536  lb_ = std::nextafter(
537  lb_, -std::numeric_limits<double>::infinity()
538  );
539  break;
540  } else {
541  lb_ = nlb;
542  adjust();
543  }
544  }
545  }
546  control_set(rhs);
547  control_check();
548  return *this;
549  }
550 
551  double inf() const {
552  return lb_;
553  }
554 
555  double sup() const {
556  return ub_;
557  }
558 
559  double estimate() const {
560  return 0.5*(lb_ + ub_);
561  }
562 
563  bool is_nan() const {
564  return !(lb_==lb_) || !(ub_==ub_);
565  }
566 
567  Sign2 sign() const {
568  if(is_nan()) {
569  return SIGN2_ERROR;
570  }
571  // Branchless (not sure it is super though...)
572  int lz = int(lb_ == 0);
573  int ln = int(lb_ < 0);
574  int lp = int(lb_ > 0);
575  int uz = int(ub_ == 0);
576  int un = int(ub_ < 0);
577  int up = int(ub_ > 0);
578  Sign2 result = Sign2(
579  ln*up*SIGN2_NP+
580  lp*up*SIGN2_PP+
581  lz*up*SIGN2_ZP+
582  ln*un*SIGN2_NN+
583  ln*uz*SIGN2_NZ
584  );
585  result = Sign2(
586  int(result) +
587  int(result==SIGN2_ZERO && !(lz&&uz)) * SIGN2_ERROR
588  );
589  return result;
590  }
591 
592  intervalRN& negate() {
593  lb_ = -lb_;
594  ub_ = -ub_;
595  std::swap(lb_, ub_);
596  control_negate();
597  control_check();
598  return *this;
599  }
600 
601  intervalRN& operator+=(const intervalRN &x) {
602  lb_ += x.lb_;
603  ub_ += x.ub_;
604  adjust();
605  control_add(x);
606  control_check();
607  return *this;
608  }
609 
610  intervalRN& operator-=(const intervalRN &x) {
611  lb_ -= x.ub_;
612  ub_ -= x.lb_;
613  adjust();
614  control_sub(x);
615  control_check();
616  return *this;
617  }
618 
619  intervalRN& operator*=(const intervalRN &x) {
620  if(!is_nan() && !x.is_nan()) {
621  double ll = lb_*x.lb_;
622  double lu = lb_*x.ub_;
623  double ul = ub_*x.lb_;
624  double uu = ub_*x.ub_;
625 
626  if(!(ll==ll)) ll = 0.0;
627  if(!(lu==lu)) lu = 0.0;
628  if(!(ul==ul)) ul = 0.0;
629  if(!(uu==uu)) uu = 0.0;
630 
631  if(lu<ll) std::swap(lu, ll);
632  if(ul<ll) std::swap(ul, ll);
633  if(uu<ll) std::swap(uu, ll);
634 
635  if(lu>uu) uu = lu;
636  if(ul>uu) uu = ul;
637 
638  lb_ = ll;
639  ub_ = uu;
640 
641  adjust();
642  } else {
643  lb_ = std::numeric_limits<double>::quiet_NaN();
644  ub_ = std::numeric_limits<double>::quiet_NaN();
645  }
646  control_mul(x);
647  control_check();
648  return *this;
649  }
650 
651  protected:
652 
653  void adjust() {
654  static constexpr double i = std::numeric_limits<double>::infinity();
655  static constexpr double e = std::numeric_limits<double>::epsilon();
656  // nextafter(1.0) - 1.0
657  static constexpr double m = std::numeric_limits<double>::min();
658  // smallest normalized
659  static constexpr double l = 1.0-e;
660  static constexpr double u = 1.0+e;
661  static constexpr double em = e*m;
662 
663  if(lb_==lb_ && ub_==ub_ && (lb_!=ub_ || (lb_!=i && lb_!=-i))) {
664 
665  if(lb_>ub_) {
666  std::swap(lb_, ub_);
667  }
668 
669  if(lb_>m) {
670  lb_ *= l;
671  } else if(lb_<-m) {
672  lb_ *= u;
673  } else {
674  lb_ -= em;
675  }
676 
677  if(ub_>m) {
678  ub_ *= u;
679  } else if(ub_<-m) {
680  ub_ *= l;
681  } else {
682  ub_ += em;
683  }
684  } else {
685  lb_ = std::numeric_limits<double>::quiet_NaN();
686  ub_ = std::numeric_limits<double>::quiet_NaN();
687  }
688  }
689 
690 #ifdef INTERVAL_CHECK
691  void control_check() {
692  intervalBase::control_check(inf(),sup());
693  }
694 #else
695  void control_check() {
696  }
697 #endif
698 
699  private:
700  double lb_;
701  double ub_;
702  };
703 
704  inline intervalRN operator+(const intervalRN& a, const intervalRN& b) {
705  intervalRN result = a;
706  return result += b;
707  }
708 
709  inline intervalRN operator-(const intervalRN& a, const intervalRN& b) {
710  intervalRN result = a;
711  return result -= b;
712  }
713 
714  inline intervalRN operator*(const intervalRN& a, const intervalRN& b) {
715  intervalRN result = a;
716  return result *= b;
717  }
718 
719  typedef intervalRN interval_nt; // Seems that valgrind does not support RU
720  //typedef intervalRU interval_nt;
721 
722 
730  template <> inline vec2Hg<interval_nt> operator-(
731  const vec2Hg<interval_nt>& p1, const vec2Hg<interval_nt>& p2
732  ) {
733  return vec2Hg<interval_nt>(
734  det2x2(p1.x,p1.w,p2.x,p2.w),
735  det2x2(p1.y,p1.w,p2.y,p2.w),
736  p1.w*p2.w
737  );
738  }
739 
740 
748  template <> inline vec3Hg<interval_nt> operator-(
749  const vec3Hg<interval_nt>& p1, const vec3Hg<interval_nt>& p2
750  ) {
751  return vec3Hg<interval_nt>(
752  det2x2(p1.x,p1.w,p2.x,p2.w),
753  det2x2(p1.y,p1.w,p2.y,p2.w),
754  det2x2(p1.z,p1.w,p2.z,p2.w),
755  p1.w * p2.w
756  );
757  }
758 
759 }
760 
761 #endif
#define geo_assert_not_reached
Sets a non reachable point in the program.
Definition: assert.h:177
#define geo_assert(x)
Verifies that a condition is met.
Definition: assert.h:149
Expansion_nt (expansion Number Type) is used to compute the sign of polynoms exactly.
Definition: expansion_nt.h:70
Base class for interval arithmetics.
Definition: interval_nt.h:71
Number type for interval arithmetics.
Definition: interval_nt.h:459
Interval arithmetics in round to upper (RU) mode.
Definition: interval_nt.h:212
2d vector with homogeneous coordinates
Definition: vechg.h:60
3d vector with homogeneous coordinates
Definition: vechg.h:185
High-level interface to multi-precision arithmetics.
Global Vorpaline namespace.
Definition: basic.h:55
Quaternion operator-(const Quaternion &a, const Quaternion &b)
Computes the difference between two Quaternion.
Definition: quaternion.h:252
void geo_argused(const T &)
Suppresses compiler warnings about unused parameters.
Definition: argused.h:60
geo_index_t index_t
The type for storing and manipulating indices.
Definition: numeric.h:329
Sign
Integer constants that represent the sign of a value.
Definition: numeric.h:68
@ ZERO
Definition: numeric.h:72
@ NEGATIVE
Definition: numeric.h:70
@ POSITIVE
Definition: numeric.h:74
T det2x2(const T &a11, const T &a12, const T &a21, const T &a22)
Computes a two-by-two determinant.
Definition: determinant.h:58
Quaternion operator+(const Quaternion &a, const Quaternion &b)
Computes the sum of two Quaternion.
Definition: quaternion.h:239
Sets FPU rounding mode for using this type of interval.
Definition: interval_nt.h:227
Generic implementation of geometric vectors in homogeneous coordinates.