40 #ifndef GEO_INTERVAL_NT
41 #define GEO_INTERVAL_NT
50 #include <xmmintrin.h>
74 static void set_FPU_round_to_nearest() {
76 _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST);
78 fesetround(FE_TONEAREST);
82 static void set_FPU_round_to_upper() {
84 _MM_SET_ROUNDING_MODE(_MM_ROUND_UP);
86 fesetround(FE_UPWARD);
100 static bool sign_is_determined(Sign2 s) {
107 static bool sign_is_non_zero(Sign2 s) {
113 static Sign convert_sign(Sign2 s) {
137 #ifdef INTERVAL_CHECK
142 control_ = x.control_;
144 void control_set(
double x) {
147 void control_negate() {
148 control_.rep().negate();
151 control_ += x.control_;
154 control_ -= x.control_;
157 control_ *= x.control_;
159 void control_check(
double inf,
double sup) {
160 typedef std::numeric_limits< double > dbl;
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;
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()
179 void control_set(
double x) {
188 void control_negate() {
199 void control_check(
double inf,
double sup) {
229 set_FPU_round_to_upper();
232 set_FPU_round_to_nearest();
252 intervalRU(
double l,
double u) : ln_(-l), u_(u) {
256 intervalRU(
const intervalRU& rhs) =
default;
258 intervalRU(
const expansion_nt& rhs) {
262 intervalRU& operator=(
const intervalRU& rhs) =
default;
264 intervalRU& operator=(
double rhs) {
272 intervalRU& operator=(
const expansion_nt& rhs) {
281 ln_ = -rhs.component(l-1);
282 u_ = rhs.component(l-1);
284 for(
int comp_idx=
int(l)-2; comp_idx>=0; --comp_idx) {
285 double comp = rhs.component(
index_t(comp_idx));
289 double new_u = u_ + comp;
292 u_, std::numeric_limits<double>::infinity()
306 double new_ln = ln_ - comp;
308 ln_ = std::nextafter(
309 ln_, std::numeric_limits<double>::infinity()
331 double estimate()
const {
333 return 0.5*(-ln_+u_);
336 bool is_nan()
const {
337 return !(ln_==ln_) || !(u_==u_);
342 int lz = int(ln_ == 0);
343 int ln = int(ln_ > 0);
344 int lp = int(ln_ < 0);
345 int uz = int(u_ == 0);
346 int un = int(u_ < 0);
347 int up = int(u_ > 0);
348 Sign2 result = Sign2(
357 int(result==SIGN2_ZERO && !(lz&&uz)) * SIGN2_ERROR
362 intervalRU& negate() {
369 intervalRU& operator+=(
const intervalRU &x) {
378 intervalRU& operator-=(
const intervalRU &x) {
388 intervalRU& operator*=(
const intervalRU &b) {
397 double lln = (-aln)*bln;
400 double uun = (-au)*bu;
404 double lu = (-aln)*bu;
405 double ul = au*(-bln);
408 ln_ = std::max(std::max(lln,lun),std::max(uln,uun));
409 u_ = std::max(std::max(ll,lu),std::max(ul,uu));
418 #ifdef INTERVAL_CHECK
419 void control_check() {
422 set_FPU_round_to_nearest();
423 intervalBase::control_check(inf(),sup());
424 set_FPU_round_to_upper();
427 void control_check() {
435 inline intervalRU
operator+(
const intervalRU& a,
const intervalRU& b) {
436 intervalRU result = a;
440 inline intervalRU
operator-(
const intervalRU& a,
const intervalRU& b) {
441 intervalRU result = a;
445 inline intervalRU operator*(
const intervalRU& a,
const intervalRU& b) {
446 intervalRU result = a;
487 intervalRN(
double l,
double u) : lb_(l), ub_(u) {
492 intervalRN(
const intervalRN& rhs) =
default;
494 intervalRN(
const expansion_nt& rhs) {
498 intervalRN& operator=(
const intervalRN& rhs) =
default;
500 intervalRN& operator=(
double rhs) {
508 intervalRN& operator=(
const expansion_nt& rhs) {
517 lb_ = rhs.component(l-1);
518 ub_ = rhs.component(l-1);
520 for(
int comp_idx=
int(l)-2; comp_idx>=0; --comp_idx) {
521 double comp = rhs.component(
index_t(comp_idx));
523 double nub = ub_ + comp;
525 ub_ = std::nextafter(
526 ub_, std::numeric_limits<double>::infinity()
534 double nlb = lb_ + comp;
536 lb_ = std::nextafter(
537 lb_, -std::numeric_limits<double>::infinity()
559 double estimate()
const {
560 return 0.5*(lb_ + ub_);
563 bool is_nan()
const {
564 return !(lb_==lb_) || !(ub_==ub_);
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(
587 int(result==SIGN2_ZERO && !(lz&&uz)) * SIGN2_ERROR
592 intervalRN& negate() {
601 intervalRN& operator+=(
const intervalRN &x) {
610 intervalRN& operator-=(
const intervalRN &x) {
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_;
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;
631 if(lu<ll) std::swap(lu, ll);
632 if(ul<ll) std::swap(ul, ll);
633 if(uu<ll) std::swap(uu, ll);
643 lb_ = std::numeric_limits<double>::quiet_NaN();
644 ub_ = std::numeric_limits<double>::quiet_NaN();
654 static constexpr
double i = std::numeric_limits<double>::infinity();
655 static constexpr
double e = std::numeric_limits<double>::epsilon();
657 static constexpr
double m = std::numeric_limits<double>::min();
659 static constexpr
double l = 1.0-e;
660 static constexpr
double u = 1.0+e;
661 static constexpr
double em = e*m;
663 if(lb_==lb_ && ub_==ub_ && (lb_!=ub_ || (lb_!=i && lb_!=-i))) {
685 lb_ = std::numeric_limits<double>::quiet_NaN();
686 ub_ = std::numeric_limits<double>::quiet_NaN();
690 #ifdef INTERVAL_CHECK
691 void control_check() {
692 intervalBase::control_check(inf(),sup());
695 void control_check() {
704 inline intervalRN
operator+(
const intervalRN& a,
const intervalRN& b) {
705 intervalRN result = a;
709 inline intervalRN
operator-(
const intervalRN& a,
const intervalRN& b) {
710 intervalRN result = a;
714 inline intervalRN operator*(
const intervalRN& a,
const intervalRN& b) {
715 intervalRN result = a;
719 typedef intervalRN interval_nt;
734 det2x2(p1.x,p1.w,p2.x,p2.w),
735 det2x2(p1.y,p1.w,p2.y,p2.w),
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),
#define geo_assert_not_reached
Sets a non reachable point in the program.
#define geo_assert(x)
Verifies that a condition is met.
Expansion_nt (expansion Number Type) is used to compute the sign of polynoms exactly.
Base class for interval arithmetics.
Number type for interval arithmetics.
Interval arithmetics in round to upper (RU) mode.
2d vector with homogeneous coordinates
3d vector with homogeneous coordinates
High-level interface to multi-precision arithmetics.
Global Vorpaline namespace.
Quaternion operator-(const Quaternion &a, const Quaternion &b)
Computes the difference between two Quaternion.
void geo_argused(const T &)
Suppresses compiler warnings about unused parameters.
geo_index_t index_t
The type for storing and manipulating indices.
Sign
Integer constants that represent the sign of a value.
T det2x2(const T &a11, const T &a12, const T &a21, const T &a22)
Computes a two-by-two determinant.
Quaternion operator+(const Quaternion &a, const Quaternion &b)
Computes the sum of two Quaternion.
Sets FPU rounding mode for using this type of interval.
Generic implementation of geometric vectors in homogeneous coordinates.