Scippy

SCIP

Solving Constraint Integer Programs

intervalarith.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file intervalarith.c
17  * @ingroup OTHER_CFILES
18  * @brief interval arithmetics for provable bounds
19  * @author Tobias Achterberg
20  * @author Stefan Vigerske
21  * @author Kati Wolter
22  */
23 
24 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
25 
26 #define _USE_MATH_DEFINES /* to get M_PI on Windows */ /*lint !750 */
27 
28 #include <stdlib.h>
29 #include <assert.h>
30 #include <math.h>
31 
32 #include "scip/def.h"
33 #include "scip/intervalarith.h"
34 #include "scip/pub_message.h"
35 #include "scip/misc.h"
36 
37 
38 /* some static initializations that need to come before enabling fenv_access
39  * (MSVC doesn't consider something like 1.5*M_PI a constant initialization if fenv_access is enabled)
40  */
41 /* first one is 1 so even indices are the maximum points */
42 static SCIP_Real sin_extremepoints[] = {M_PI_2, 1.5*M_PI, 2.5*M_PI, 3.5*M_PI};
43 /* first one is -1 so even indices are the minimum points */
44 static SCIP_Real cos_extremepoints[] = {M_PI, 2*M_PI, 3*M_PI};
45 
46 /* Inform compiler that this code accesses the floating-point environment, so that
47  * certain optimizations should be omitted (http://www.cplusplus.com/reference/cfenv/FENV_ACCESS/).
48  * Not supported by Clang (gives warning) and GCC (silently), at the moment.
49  */
50 #if defined(__INTEL_COMPILER) || defined(_MSC_VER)
51 #pragma fenv_access (on)
52 #elif defined __GNUC__
53 #pragma STDC FENV_ACCESS ON
54 #endif
55 
56 /* Unfortunately, the FENV_ACCESS pragma is essentially ignored by GCC at the moment (2019),
57  * see #2650 and https://gcc.gnu.org/bugzilla/show_bug.cgi?id=34678.
58  * There are ways to work around this by declaring variables volatile or inserting more assembler code,
59  * but there is always the danger that something would be overlooked.
60  * A more drastic but safer way seems to be to just disable all compiler optimizations for this file.
61  * The Intel compiler seems to implement FENV_ACCESS correctly, but also defines __GNUC__.
62  */
63 #if defined(__GNUC__) && !defined( __INTEL_COMPILER)
64 #pragma GCC push_options
65 #pragma GCC optimize ("O0")
66 #endif
67 
68 #ifdef SCIP_ROUNDING_FE
69 #define ROUNDING
70 /*
71  * Linux rounding operations
72  */
73 
74 #include <fenv.h>
75 
76 /** Linux rounding mode settings */
77 #define SCIP_ROUND_DOWNWARDS FE_DOWNWARD /**< round always down */
78 #define SCIP_ROUND_UPWARDS FE_UPWARD /**< round always up */
79 #define SCIP_ROUND_NEAREST FE_TONEAREST /**< round always to nearest */
80 #define SCIP_ROUND_ZERO FE_TOWARDZERO /**< round always towards 0.0 */
81 
82 /** returns whether rounding mode control is available */
84  void
85  )
86 {
87  return TRUE;
88 }
89 
90 /** sets rounding mode of floating point operations */
91 static
93  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
94  )
95 {
96 #ifndef NDEBUG
97  if( fesetround(roundmode) != 0 )
98  {
99  SCIPerrorMessage("error setting rounding mode to %d\n", roundmode);
100  abort();
101  }
102 #else
103  (void) fesetround(roundmode);
104 #endif
105 }
106 
107 /** gets current rounding mode of floating point operations */
108 static
110  void
111  )
112 {
113  return (SCIP_ROUNDMODE)fegetround();
114 }
115 
116 #endif
117 
118 
119 
120 #ifdef SCIP_ROUNDING_FP
121 #define ROUNDING
122 /*
123  * OSF rounding operations
124  */
125 
126 #include <float.h>
127 
128 /** OSF rounding mode settings */
129 #define SCIP_ROUND_DOWNWARDS FP_RND_RM /**< round always down */
130 #define SCIP_ROUND_UPWARDS FP_RND_RP /**< round always up */
131 #define SCIP_ROUND_NEAREST FP_RND_RN /**< round always to nearest */
132 #define SCIP_ROUND_ZERO FP_RND_RZ /**< round always towards 0.0 */
133 
134 /** returns whether rounding mode control is available */
136  void
137  )
138 {
139  return TRUE;
140 }
141 
142 /** sets rounding mode of floating point operations */
143 static
145  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
146  )
147 {
148 #ifndef NDEBUG
149  if( write_rnd(roundmode) != 0 )
150  {
151  SCIPerrorMessage("error setting rounding mode to %d\n", roundmode);
152  abort();
153  }
154 #else
155  (void) write_rnd(roundmode);
156 #endif
157 }
158 
159 /** gets current rounding mode of floating point operations */
160 static
162  void
163  )
164 {
165  return read_rnd();
166 }
167 
168 #endif
169 
170 
171 
172 #ifdef SCIP_ROUNDING_MS
173 #define ROUNDING
174 /*
175  * Microsoft compiler rounding operations
176  */
177 
178 #include <float.h>
179 
180 /** Microsoft rounding mode settings */
181 #define SCIP_ROUND_DOWNWARDS RC_DOWN /**< round always down */
182 #define SCIP_ROUND_UPWARDS RC_UP /**< round always up */
183 #define SCIP_ROUND_NEAREST RC_NEAR /**< round always to nearest */
184 #define SCIP_ROUND_ZERO RC_CHOP /**< round always towards zero */
185 
186 /** returns whether rounding mode control is available */
188  void
189  )
190 {
191  return TRUE;
192 }
193 
194 /** sets rounding mode of floating point operations */
195 static
197  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
198  )
199 {
200 #ifndef NDEBUG
201  if( (_controlfp(roundmode, _MCW_RC) & _MCW_RC) != roundmode )
202  {
203  SCIPerrorMessage("error setting rounding mode to %x\n", roundmode);
204  abort();
205  }
206 #else
207  (void) _controlfp(roundmode, _MCW_RC);
208 #endif
209 }
210 
211 /** gets current rounding mode of floating point operations */
212 static
214  void
215  )
216 {
217  return _controlfp(0, 0) & _MCW_RC;
218 }
219 #endif
220 
221 
222 
223 #ifndef ROUNDING
224 /*
225  * rouding operations not available
226  */
227 #define SCIP_ROUND_DOWNWARDS 0 /**< round always down */
228 #define SCIP_ROUND_UPWARDS 1 /**< round always up */
229 #define SCIP_ROUND_NEAREST 2 /**< round always to nearest */
230 #define SCIP_ROUND_ZERO 3 /**< round always towards zero */
231 
232 /** returns whether rounding mode control is available */
234  void
235  )
236 {
237  return FALSE;
238 }
239 
240 /** sets rounding mode of floating point operations */ /*lint -e715*/
241 static
243  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
244  )
245 { /*lint --e{715}*/
246  SCIPerrorMessage("setting rounding mode not available - interval arithmetic is invalid!\n");
247 }
248 
249 /** gets current rounding mode of floating point operations */
250 static
252  void
253  )
254 {
255  return SCIP_ROUND_NEAREST;
256 }
257 #else
258 #undef ROUNDING
259 #endif
260 
261 /** sets rounding mode of floating point operations */
263  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
264  )
265 {
266  intervalSetRoundingMode(roundmode);
267 }
268 
269 /** gets current rounding mode of floating point operations */
271  void
272  )
273 {
274  return intervalGetRoundingMode();
275 }
276 
277 #if defined(__GNUC__) && (defined(__i386__) || defined(__x86_64__)) /* gcc or icc compiler on x86 32bit or 64bit */
278 
279 /** gets the negation of a double
280  * Do this in a way that the compiler does not "optimize" it away, which usually does not considers rounding modes.
281  * However, compiling with -frounding-math would allow to return -x here.
282  * @todo We now set the FENV_ACCESS pragma to on, which is the same as -frounding-math, so we might be able to eliminate this.
283  */
284 static
285 double negate(
286  /* we explicitly use double here, since I'm not sure the assembler code would work as it for other float's */
287  double x /**< number that should be negated */
288  )
289 {
290  /* The following line of code is taken from GAOL, http://sourceforge.net/projects/gaol. */
291  __asm volatile ("fldl %1; fchs; fstpl %0" : "=m" (x) : "m" (x));
292  return x;
293 }
294 
295 /* cl or icl compiler on 32bit windows or icl compiler on 64bit windows
296  * cl on 64bit windows does not seem to support inline assembler
297  */
298 #elif defined(_MSC_VER) && (defined(__INTEL_COMPILER) || !defined(_M_X64))
299 
300 /** gets the negation of a double
301  * Do this in a way that the compiler does not "optimize" it away, which usually does not considers rounding modes.
302  */
303 static
304 double negate(
305  /* we explicitly use double here, since I'm not sure the assembler code would work as it for other float's */
306  double x /**< number that should be negated */
307  )
308 {
309  /* The following lines of code are taken from GAOL, http://sourceforge.net/projects/gaol. */
310  __asm {
311  fld x
312  fchs
313  fstp x
314  }
315  return x;
316 }
317 
318 #else /* unknown compiler or MSVS 64bit */
319 
320 /** gets the negation of a double
321  *
322  * Fallback implementation that calls the negation method from misc.o.
323  * Having the implementation in a different object file will hopefully prevent
324  * it from being "optimized away".
325  */
326 static
328  SCIP_Real x /**< number that should be negated */
329  )
330 {
331  return SCIPnegateReal(x);
332 }
333 
334 #endif
335 
336 
337 /** sets rounding mode of floating point operations to downwards rounding */
339  void
340  )
341 {
343 }
344 
345 /** sets rounding mode of floating point operations to upwards rounding */
347  void
348  )
349 {
351 }
352 
353 /** sets rounding mode of floating point operations to nearest rounding */
355  void
356  )
357 {
359 }
360 
361 /** sets rounding mode of floating point operations to towards zero rounding */
363  void
364  )
365 {
367 }
368 
369 /** negates a number in a way that the compiler does not optimize it away */
371  SCIP_Real x /**< number to negate */
372  )
373 {
374  return negate((double)x);
375 }
376 
377 /*
378  * Interval arithmetic operations
379  */
380 
381 /* In debug mode, the following methods are implemented as function calls to ensure
382  * type validity.
383  * In optimized mode, the methods are implemented as defines to improve performance.
384  * However, we want to have them in the library anyways, so we have to undef the defines.
385  */
386 
387 #undef SCIPintervalGetInf
388 #undef SCIPintervalGetSup
389 #undef SCIPintervalSet
390 #undef SCIPintervalSetBounds
391 #undef SCIPintervalSetEmpty
392 #undef SCIPintervalIsEmpty
393 #undef SCIPintervalSetEntire
394 #undef SCIPintervalIsEntire
395 #undef SCIPintervalIsPositiveInfinity
396 #undef SCIPintervalIsNegativeInfinity
397 
398 /** returns infimum of interval */
400  SCIP_INTERVAL interval /**< interval */
401  )
402 {
403  return interval.inf;
404 }
405 
406 /** returns supremum of interval */
408  SCIP_INTERVAL interval /**< interval */
409  )
410 {
411  return interval.sup;
412 }
413 
414 /** stores given value as interval */
416  SCIP_INTERVAL* resultant, /**< interval to store value into */
417  SCIP_Real value /**< value to store */
418  )
419 {
420  assert(resultant != NULL);
421 
422  resultant->inf = value;
423  resultant->sup = value;
424 }
425 
426 /** stores given infimum and supremum as interval */
428  SCIP_INTERVAL* resultant, /**< interval to store value into */
429  SCIP_Real inf, /**< value to store as infimum */
430  SCIP_Real sup /**< value to store as supremum */
431  )
432 {
433  assert(resultant != NULL);
434  assert(inf <= sup);
435 
436  resultant->inf = inf;
437  resultant->sup = sup;
438 }
439 
440 /** sets interval to empty interval, which will be [infinity, -infinity] */
442  SCIP_INTERVAL* resultant /**< resultant interval of operation */
443  )
444 {
445  assert(resultant != NULL);
446 
447  resultant->inf = 1.0;
448  resultant->sup = -1.0;
449 }
450 
451 /** indicates whether interval is empty, i.e., whether inf > sup */
453  SCIP_Real infinity, /**< value for infinity */
454  SCIP_INTERVAL operand /**< operand of operation */
455  )
456 {
457  if( operand.sup >= infinity || operand.inf <= -infinity )
458  return FALSE;
459 
460  return operand.sup < operand.inf;
461 }
462 
463 /** sets interval to entire [-infinity, +infinity] */
465  SCIP_Real infinity, /**< value for infinity */
466  SCIP_INTERVAL* resultant /**< resultant interval of operation */
467  )
468 {
469  assert(resultant != NULL);
470 
471  resultant->inf = -infinity;
472  resultant->sup = infinity;
473 }
474 
475 /** indicates whether interval is entire, i.e., whether inf <= -infinity and sup >= infinity */
477  SCIP_Real infinity, /**< value for infinity */
478  SCIP_INTERVAL operand /**< operand of operation */
479  )
480 {
481  return operand.inf <= -infinity && operand.sup >= infinity;
482 }
483 
484 /** indicates whether interval is positive infinity, i.e., [infinity, infinity] */
486  SCIP_Real infinity, /**< value for infinity */
487  SCIP_INTERVAL operand /**< operand of operation */
488  )
489 {
490  return operand.inf >= infinity && operand.sup >= operand.inf;
491 }
492 
493 /** indicates whether interval is negative infinity, i.e., [-infinity, -infinity] */
495  SCIP_Real infinity, /**< value for infinity */
496  SCIP_INTERVAL operand /**< operand of operation */
497  )
498 {
499  return operand.sup <= -infinity && operand.inf <= operand.sup;
500 }
501 
502 /** indicates whether operand1 is contained in operand2 */
504  SCIP_Real infinity, /**< value for infinity */
505  SCIP_INTERVAL operand1, /**< first operand of operation */
506  SCIP_INTERVAL operand2 /**< second operand of operation */
507  )
508 {
509  /* the empty interval is contained everywhere */
510  if( operand1.inf > operand1.sup )
511  return TRUE;
512 
513  /* something not-empty is not contained in the empty interval */
514  if( operand2.inf > operand2.sup )
515  return FALSE;
516 
517  return (MAX(-infinity, operand1.inf) >= operand2.inf) &&
518  ( MIN( infinity, operand1.sup) <= operand2.sup);
519 }
520 
521 /** indicates whether operand1 and operand2 are disjoint */
523  SCIP_INTERVAL operand1, /**< first operand of operation */
524  SCIP_INTERVAL operand2 /**< second operand of operation */
525  )
526 {
527  return (operand1.sup < operand2.inf) || (operand2.sup < operand1.inf);
528 }
529 
530 /** intersection of two intervals */
532  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
533  SCIP_INTERVAL operand1, /**< first operand of operation */
534  SCIP_INTERVAL operand2 /**< second operand of operation */
535  )
536 {
537  assert(resultant != NULL);
538 
539  resultant->inf = MAX(operand1.inf, operand2.inf); /*lint !e644*/
540  resultant->sup = MIN(operand1.sup, operand2.sup); /*lint !e644*/
541 }
542 
543 /** interval enclosure of the union of two intervals */
545  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
546  SCIP_INTERVAL operand1, /**< first operand of operation */
547  SCIP_INTERVAL operand2 /**< second operand of operation */
548  )
549 {
550  assert(resultant != NULL);
551 
552  if( operand1.inf > operand1.sup )
553  {
554  /* operand1 is empty */
555  *resultant = operand2;
556  return;
557  }
558 
559  if( operand2.inf > operand2.sup )
560  {
561  /* operand2 is empty */
562  *resultant = operand1;
563  return;
564  }
565 
566  resultant->inf = MIN(operand1.inf, operand2.inf);
567  resultant->sup = MAX(operand1.sup, operand2.sup);
568 }
569 
570 /** adds operand1 and operand2 and stores infimum of result in infimum of resultant */
572  SCIP_Real infinity, /**< value for infinity */
573  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
574  SCIP_INTERVAL operand1, /**< first operand of operation */
575  SCIP_INTERVAL operand2 /**< second operand of operation */
576  )
577 {
579  assert(resultant != NULL);
580 
581  /* [a,...] + [-inf,...] = [-inf,...] for all a, in particular, [+inf,...] + [-inf,...] = [-inf,...] */
582  if( operand1.inf <= -infinity || operand2.inf <= -infinity )
583  {
584  resultant->inf = -infinity;
585  }
586  /* [a,...] + [+inf,...] = [+inf,...] for all a > -inf */
587  else if( operand1.inf >= infinity || operand2.inf >= infinity )
588  {
589  resultant->inf = infinity;
590  }
591  else
592  {
593  resultant->inf = operand1.inf + operand2.inf;
594  }
595 }
596 
597 /** adds operand1 and operand2 and stores supremum of result in supremum of resultant */
599  SCIP_Real infinity, /**< value for infinity */
600  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
601  SCIP_INTERVAL operand1, /**< first operand of operation */
602  SCIP_INTERVAL operand2 /**< second operand of operation */
603  )
604 {
606  assert(resultant != NULL);
607 
608  /* [...,b] + [...,+inf] = [...,+inf] for all b, in particular, [...,-inf] + [...,+inf] = [...,+inf] */
609  if( operand1.sup >= infinity || operand2.sup >= infinity )
610  {
611  resultant->sup = infinity;
612  }
613  /* [...,b] + [...,-inf] = [...,-inf] for all b < +inf */
614  else if( operand1.sup <= -infinity || operand2.sup <= -infinity )
615  {
616  resultant->sup = -infinity;
617  }
618  else
619  {
620  resultant->sup = operand1.sup + operand2.sup;
621  }
622 }
623 
624 /** adds operand1 and operand2 and stores result in resultant */
626  SCIP_Real infinity, /**< value for infinity */
627  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
628  SCIP_INTERVAL operand1, /**< first operand of operation */
629  SCIP_INTERVAL operand2 /**< second operand of operation */
630  )
631 {
632  SCIP_ROUNDMODE roundmode;
633 
634  assert(resultant != NULL);
635  assert(!SCIPintervalIsEmpty(infinity, operand1));
636  assert(!SCIPintervalIsEmpty(infinity, operand2));
637 
638  roundmode = intervalGetRoundingMode();
639 
640  /* compute infimum of result */
642  SCIPintervalAddInf(infinity, resultant, operand1, operand2);
643 
644  /* compute supremum of result */
646  SCIPintervalAddSup(infinity, resultant, operand1, operand2);
647 
648  intervalSetRoundingMode(roundmode);
649 }
650 
651 /** adds operand1 and scalar operand2 and stores result in resultant */
653  SCIP_Real infinity, /**< value for infinity */
654  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
655  SCIP_INTERVAL operand1, /**< first operand of operation */
656  SCIP_Real operand2 /**< second operand of operation */
657  )
658 { /*lint --e{644}*/
659  SCIP_ROUNDMODE roundmode;
660 
661  assert(resultant != NULL);
662  assert(!SCIPintervalIsEmpty(infinity, operand1));
663 
664  roundmode = intervalGetRoundingMode();
665 
666  /* -inf + something >= -inf */
667  if( operand1.inf <= -infinity || operand2 <= -infinity )
668  {
669  resultant->inf = -infinity;
670  }
671  else if( operand1.inf >= infinity || operand2 >= infinity )
672  {
673  /* inf + finite = inf, inf + inf = inf */
674  resultant->inf = infinity;
675  }
676  else
677  {
679  resultant->inf = operand1.inf + operand2;
680  }
681 
682  /* inf + something <= inf */
683  if( operand1.sup >= infinity || operand2 >= infinity )
684  {
685  resultant->sup = infinity;
686  }
687  else if( operand1.sup <= -infinity || operand2 <= -infinity )
688  {
689  /* -inf + finite = -inf, -inf + (-inf) = -inf */
690  resultant->sup = -infinity;
691  }
692  else
693  {
695  resultant->sup = operand1.sup + operand2;
696  }
697 
698  intervalSetRoundingMode(roundmode);
699 }
700 
701 /** adds vector operand1 and vector operand2 and stores result in vector resultant */
703  SCIP_Real infinity, /**< value for infinity */
704  SCIP_INTERVAL* resultant, /**< array of resultant intervals of operation */
705  int length, /**< length of arrays */
706  SCIP_INTERVAL* operand1, /**< array of first operands of operation */
707  SCIP_INTERVAL* operand2 /**< array of second operands of operation */
708  )
709 {
710  SCIP_ROUNDMODE roundmode;
711  int i;
712 
713  roundmode = intervalGetRoundingMode();
714 
715  /* compute infimums of resultant array */
717  for( i = 0; i < length; ++i )
718  {
719  SCIPintervalAddInf(infinity, &resultant[i], operand1[i], operand2[i]);
720  }
721  /* compute supremums of result array */
723  for( i = 0; i < length; ++i )
724  {
725  SCIPintervalAddSup(infinity, &resultant[i], operand1[i], operand2[i]);
726  }
727 
728  intervalSetRoundingMode(roundmode);
729 }
730 
731 /** subtracts operand2 from operand1 and stores result in resultant */
733  SCIP_Real infinity, /**< value for infinity */
734  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
735  SCIP_INTERVAL operand1, /**< first operand of operation */
736  SCIP_INTERVAL operand2 /**< second operand of operation */
737  )
738 { /*lint --e{644}*/
739  SCIP_ROUNDMODE roundmode;
740 
741  assert(resultant != NULL);
742  assert(!SCIPintervalIsEmpty(infinity, operand1));
743  assert(!SCIPintervalIsEmpty(infinity, operand2));
744 
745  roundmode = intervalGetRoundingMode();
746 
747  if( operand1.inf <= -infinity || operand2.sup >= infinity )
748  resultant->inf = -infinity;
749  /* [a,b] - [-inf,-inf] = [+inf,+inf] */
750  else if( operand1.inf >= infinity || operand2.sup <= -infinity )
751  {
752  resultant->inf = infinity;
753  resultant->sup = infinity;
754  return;
755  }
756  else
757  {
759  resultant->inf = operand1.inf - operand2.sup;
760  }
761 
762  if( operand1.sup >= infinity || operand2.inf <= -infinity )
763  resultant->sup = infinity;
764  /* [a,b] - [+inf,+inf] = [-inf,-inf] */
765  else if( operand1.sup <= -infinity || operand2.inf >= infinity )
766  {
767  assert(resultant->inf == -infinity); /* should be set above, since operand1.inf <= operand1.sup <= -infinity */ /*lint !e777*/
768  resultant->sup = -infinity;
769  }
770  else
771  {
773  resultant->sup = operand1.sup - operand2.inf;
774  }
775 
776  intervalSetRoundingMode(roundmode);
777 }
778 
779 /** subtracts scalar operand2 from operand1 and stores result in resultant */
781  SCIP_Real infinity, /**< value for infinity */
782  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
783  SCIP_INTERVAL operand1, /**< first operand of operation */
784  SCIP_Real operand2 /**< second operand of operation */
785  )
786 {
787  SCIPintervalAddScalar(infinity, resultant, operand1, -operand2);
788 }
789 
790 /** multiplies operand1 with operand2 and stores infimum of result in infimum of resultant */
792  SCIP_Real infinity, /**< value for infinity */
793  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
794  SCIP_INTERVAL operand1, /**< first operand of operation; can be +/-inf */
795  SCIP_INTERVAL operand2 /**< second operand of operation; can be +/-inf */
796  )
797 {
798  assert(resultant != NULL);
799  assert(!SCIPintervalIsEmpty(infinity, operand1));
800  assert(!SCIPintervalIsEmpty(infinity, operand2));
801 
803 
804  if( operand1.inf >= infinity )
805  {
806  /* operand1 is infinity scalar */
807  assert(operand1.sup >= infinity);
808  SCIPintervalMulScalarInf(infinity, resultant, operand2, infinity);
809  }
810  else if( operand2.inf >= infinity )
811  {
812  /* operand2 is infinity scalar */
813  assert(operand2.sup >= infinity);
814  SCIPintervalMulScalarInf(infinity, resultant, operand1, infinity);
815  }
816  else if( operand1.sup <= -infinity )
817  {
818  /* operand1 is -infinity scalar */
819  assert(operand1.inf <= -infinity);
820  SCIPintervalMulScalarInf(infinity, resultant, operand2, -infinity);
821  }
822  else if( operand2.sup <= -infinity )
823  {
824  /* operand2 is -infinity scalar */
825  assert(operand2.inf <= -infinity);
826  SCIPintervalMulScalarInf(infinity, resultant, operand1, -infinity);
827  }
828  else if( ( operand1.inf <= -infinity && operand2.sup > 0.0 )
829  || ( operand1.sup > 0.0 && operand2.inf <= -infinity )
830  || ( operand1.inf < 0.0 && operand2.sup >= infinity )
831  || ( operand1.sup >= infinity && operand2.inf < 0.0 ) )
832  {
833  resultant->inf = -infinity;
834  }
835  else
836  {
837  SCIP_Real cand1;
838  SCIP_Real cand2;
839  SCIP_Real cand3;
840  SCIP_Real cand4;
841 
842  cand1 = operand1.inf * operand2.inf;
843  cand2 = operand1.inf * operand2.sup;
844  cand3 = operand1.sup * operand2.inf;
845  cand4 = operand1.sup * operand2.sup;
846  resultant->inf = MIN(MIN(cand1, cand2), MIN(cand3, cand4));
847  }
848 }
849 
850 /** multiplies operand1 with operand2 and stores supremum of result in supremum of resultant */
852  SCIP_Real infinity, /**< value for infinity */
853  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
854  SCIP_INTERVAL operand1, /**< first operand of operation; can be +/-inf */
855  SCIP_INTERVAL operand2 /**< second operand of operation; can be +/-inf */
856  )
857 {
858  assert(resultant != NULL);
859  assert(!SCIPintervalIsEmpty(infinity, operand1));
860  assert(!SCIPintervalIsEmpty(infinity, operand2));
861 
863 
864  if( operand1.inf >= infinity )
865  {
866  /* operand1 is infinity scalar */
867  assert(operand1.sup >= infinity);
868  SCIPintervalMulScalarSup(infinity, resultant, operand2, infinity);
869  }
870  else if( operand2.inf >= infinity )
871  {
872  /* operand2 is infinity scalar */
873  assert(operand2.sup >= infinity);
874  SCIPintervalMulScalarSup(infinity, resultant, operand1, infinity);
875  }
876  else if( operand1.sup <= -infinity )
877  {
878  /* operand1 is -infinity scalar */
879  assert(operand1.inf <= -infinity);
880  SCIPintervalMulScalarSup(infinity, resultant, operand2, -infinity);
881  }
882  else if( operand2.sup <= -infinity )
883  {
884  /* operand2 is -infinity scalar */
885  assert(operand2.inf <= -infinity);
886  SCIPintervalMulScalarSup(infinity, resultant, operand1, -infinity);
887  }
888  else if( ( operand1.inf <= -infinity && operand2.inf < 0.0 )
889  || ( operand1.inf < 0.0 && operand2.inf <= -infinity )
890  || ( operand1.sup > 0.0 && operand2.sup >= infinity )
891  || ( operand1.sup >= infinity && operand2.sup > 0.0 ) )
892  {
893  resultant->sup = infinity;
894  }
895  else
896  {
897  SCIP_Real cand1;
898  SCIP_Real cand2;
899  SCIP_Real cand3;
900  SCIP_Real cand4;
901 
902  cand1 = operand1.inf * operand2.inf;
903  cand2 = operand1.inf * operand2.sup;
904  cand3 = operand1.sup * operand2.inf;
905  cand4 = operand1.sup * operand2.sup;
906  resultant->sup = MAX(MAX(cand1, cand2), MAX(cand3, cand4));
907  }
908 }
909 
910 /** multiplies operand1 with operand2 and stores result in resultant */
912  SCIP_Real infinity, /**< value for infinity */
913  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
914  SCIP_INTERVAL operand1, /**< first operand of operation; can be +/-inf */
915  SCIP_INTERVAL operand2 /**< second operand of operation; can be +/-inf */
916  )
917 {
918  SCIP_ROUNDMODE roundmode;
919 
920  assert(resultant != NULL);
921  assert(!SCIPintervalIsEmpty(infinity, operand1));
922  assert(!SCIPintervalIsEmpty(infinity, operand2));
923 
924  roundmode = intervalGetRoundingMode();
925 
926  /* compute infimum result */
928  SCIPintervalMulInf(infinity, resultant, operand1, operand2);
929 
930  /* compute supremum of result */
932  SCIPintervalMulSup(infinity, resultant, operand1, operand2);
933 
934  intervalSetRoundingMode(roundmode);
935 }
936 
937 /** multiplies operand1 with scalar operand2 and stores infimum of result in infimum of resultant */
939  SCIP_Real infinity, /**< value for infinity */
940  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
941  SCIP_INTERVAL operand1, /**< first operand of operation */
942  SCIP_Real operand2 /**< second operand of operation; can be +/- inf */
943  )
944 {
946  assert(resultant != NULL);
947  assert(!SCIPintervalIsEmpty(infinity, operand1));
948 
949  if( operand2 >= infinity )
950  {
951  /* result.inf defined by sign of operand1.inf */
952  if( operand1.inf > 0 )
953  resultant->inf = infinity;
954  else if( operand1.inf < 0 )
955  resultant->inf = -infinity;
956  else
957  resultant->inf = 0.0;
958  }
959  else if( operand2 <= -infinity )
960  {
961  /* result.inf defined by sign of operand1.sup */
962  if( operand1.sup > 0 )
963  resultant->inf = -infinity;
964  else if( operand1.sup < 0 )
965  resultant->inf = infinity;
966  else
967  resultant->inf = 0.0;
968  }
969  else if( operand2 == 0.0 )
970  {
971  resultant->inf = 0.0;
972  }
973  else if( operand2 > 0.0 )
974  {
975  if( operand1.inf <= -infinity )
976  resultant->inf = -infinity;
977  else if( operand1.inf >= infinity )
978  resultant->inf = infinity;
979  else
980  resultant->inf = operand1.inf * operand2;
981  }
982  else
983  {
984  if( operand1.sup >= infinity )
985  resultant->inf = -infinity;
986  else if( operand1.sup <= -infinity )
987  resultant->inf = infinity;
988  else
989  resultant->inf = operand1.sup * operand2;
990  }
991 }
992 
993 /** multiplies operand1 with scalar operand2 and stores supremum of result in supremum of resultant */
995  SCIP_Real infinity, /**< value for infinity */
996  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
997  SCIP_INTERVAL operand1, /**< first operand of operation */
998  SCIP_Real operand2 /**< second operand of operation; can be +/- inf */
999  )
1000 {
1002  assert(resultant != NULL);
1003  assert(!SCIPintervalIsEmpty(infinity, operand1));
1004 
1005  if( operand2 >= infinity )
1006  {
1007  /* result.sup defined by sign of operand1.sup */
1008  if( operand1.sup > 0 )
1009  resultant->sup = infinity;
1010  else if( operand1.sup < 0 )
1011  resultant->sup = -infinity;
1012  else
1013  resultant->sup = 0.0;
1014  }
1015  else if( operand2 <= -infinity )
1016  {
1017  /* result.sup defined by sign of operand1.inf */
1018  if( operand1.inf > 0 )
1019  resultant->sup = -infinity;
1020  else if( operand1.inf < 0 )
1021  resultant->sup = infinity;
1022  else
1023  resultant->sup = 0.0;
1024  }
1025  else if( operand2 == 0.0 )
1026  {
1027  resultant->sup = 0.0;
1028  }
1029  else if( operand2 > 0.0 )
1030  {
1031  if( operand1.sup >= infinity )
1032  resultant->sup = infinity;
1033  else if( operand1.sup <= -infinity )
1034  resultant->sup = -infinity;
1035  else
1036  resultant->sup = operand1.sup * operand2;
1037  }
1038  else
1039  {
1040  if( operand1.inf <= -infinity )
1041  resultant->sup = infinity;
1042  else if( operand1.inf >= infinity )
1043  resultant->sup = -infinity;
1044  else
1045  resultant->sup = operand1.inf * operand2;
1046  }
1047 }
1048 
1049 /** multiplies operand1 with scalar operand2 and stores result in resultant */
1051  SCIP_Real infinity, /**< value for infinity */
1052  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1053  SCIP_INTERVAL operand1, /**< first operand of operation */
1054  SCIP_Real operand2 /**< second operand of operation */
1055  )
1056 {
1057  SCIP_ROUNDMODE roundmode;
1058 
1059  assert(resultant != NULL);
1060  assert(!SCIPintervalIsEmpty(infinity, operand1));
1061 
1062  if( operand2 == 1.0 ) /*lint !e777*/
1063  {
1064  *resultant = operand1;
1065  return;
1066  }
1067 
1068  if( operand2 == -1.0 ) /*lint !e777*/
1069  {
1070  resultant->inf = -operand1.sup;
1071  resultant->sup = -operand1.inf;
1072  return;
1073  }
1074 
1075  roundmode = intervalGetRoundingMode();
1076 
1077  /* compute infimum result */
1079  SCIPintervalMulScalarInf(infinity, resultant, operand1, operand2);
1080 
1081  /* compute supremum of result */
1083  SCIPintervalMulScalarSup(infinity, resultant, operand1, operand2);
1084 
1085  intervalSetRoundingMode(roundmode);
1086 }
1087 
1088 /** divides operand1 by operand2 and stores result in resultant */
1090  SCIP_Real infinity, /**< value for infinity */
1091  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1092  SCIP_INTERVAL operand1, /**< first operand of operation */
1093  SCIP_INTERVAL operand2 /**< second operand of operation */
1094  )
1095 {
1096  SCIP_ROUNDMODE roundmode;
1097  SCIP_INTERVAL intmed;
1098 
1099  assert(resultant != NULL);
1100  assert(!SCIPintervalIsEmpty(infinity, operand1));
1101  assert(!SCIPintervalIsEmpty(infinity, operand2));
1102 
1103  if( operand2.inf <= 0.0 && operand2.sup >= 0.0 )
1104  { /* division by [0,0] or interval containing 0 gives [-inf, +inf] */
1105  resultant->inf = -infinity;
1106  resultant->sup = infinity;
1107  return;
1108  }
1109 
1110  if( operand1.inf == 0.0 && operand1.sup == 0.0 )
1111  { /* division of [0,0] by something nonzero */
1112  SCIPintervalSet(resultant, 0.0);
1113  return;
1114  }
1115 
1116  roundmode = intervalGetRoundingMode();
1117 
1118  /* division by nonzero: resultant = x * (1/y) */
1119  if( operand2.sup >= infinity || operand2.sup <= -infinity )
1120  {
1121  intmed.inf = 0.0;
1122  }
1123  else
1124  {
1126  intmed.inf = 1.0 / operand2.sup;
1127  }
1128  if( operand2.inf <= -infinity || operand2.inf >= infinity )
1129  {
1130  intmed.sup = 0.0;
1131  }
1132  else
1133  {
1135  intmed.sup = 1.0 / operand2.inf;
1136  }
1137  SCIPintervalMul(infinity, resultant, operand1, intmed);
1138 
1139  intervalSetRoundingMode(roundmode);
1140 }
1141 
1142 /** divides operand1 by scalar operand2 and stores result in resultant */
1144  SCIP_Real infinity, /**< value for infinity */
1145  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1146  SCIP_INTERVAL operand1, /**< first operand of operation */
1147  SCIP_Real operand2 /**< second operand of operation */
1148  )
1149 {
1150  SCIP_ROUNDMODE roundmode;
1151 
1152  assert(resultant != NULL);
1153  assert(!SCIPintervalIsEmpty(infinity, operand1));
1154 
1155  roundmode = intervalGetRoundingMode();
1156 
1157  if( operand2 >= infinity || operand2 <= -infinity )
1158  {
1159  /* division by +/-infinity is 0.0 */
1160  resultant->inf = 0.0;
1161  resultant->sup = 0.0;
1162  }
1163  else if( operand2 > 0.0 )
1164  {
1165  if( operand1.inf <= -infinity )
1166  resultant->inf = -infinity;
1167  else if( operand1.inf >= infinity )
1168  {
1169  /* infinity / + = infinity */
1170  resultant->inf = infinity;
1171  }
1172  else
1173  {
1175  resultant->inf = operand1.inf / operand2;
1176  }
1177 
1178  if( operand1.sup >= infinity )
1179  resultant->sup = infinity;
1180  else if( operand1.sup <= -infinity )
1181  {
1182  /* -infinity / + = -infinity */
1183  resultant->sup = -infinity;
1184  }
1185  else
1186  {
1188  resultant->sup = operand1.sup / operand2;
1189  }
1190  }
1191  else if( operand2 < 0.0 )
1192  {
1193  if( operand1.sup >= infinity )
1194  resultant->inf = -infinity;
1195  else if( operand1.sup <= -infinity )
1196  {
1197  /* -infinity / - = infinity */
1198  resultant->inf = infinity;
1199  }
1200  else
1201  {
1203  resultant->inf = operand1.sup / operand2;
1204  }
1205 
1206  if( operand1.inf <= -infinity )
1207  resultant->sup = infinity;
1208  else if( operand1.inf >= infinity )
1209  {
1210  /* infinity / - = -infinity */
1211  resultant->sup = -infinity;
1212  }
1213  else
1214  {
1216  resultant->sup = operand1.inf / operand2;
1217  }
1218  }
1219  else
1220  { /* division by 0.0 */
1221  if( operand1.inf >= 0 )
1222  {
1223  /* [+,+] / [0,0] = [+inf, +inf] */
1224  resultant->inf = infinity;
1225  resultant->sup = infinity;
1226  }
1227  else if( operand1.sup <= 0 )
1228  {
1229  /* [-,-] / [0,0] = [-inf, -inf] */
1230  resultant->inf = -infinity;
1231  resultant->sup = -infinity;
1232  }
1233  else
1234  {
1235  /* [-,+] / [0,0] = [-inf, +inf] */
1236  resultant->inf = -infinity;
1237  resultant->sup = infinity;
1238  }
1239  return;
1240  }
1241 
1242  intervalSetRoundingMode(roundmode);
1243 }
1244 
1245 /** computes the scalar product of two vectors of intervals and stores result in resultant */
1247  SCIP_Real infinity, /**< value for infinity */
1248  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1249  int length, /**< length of vectors */
1250  SCIP_INTERVAL* operand1, /**< first vector as array of intervals; can have +/-inf entries */
1251  SCIP_INTERVAL* operand2 /**< second vector as array of intervals; can have +/-inf entries */
1252  )
1253 {
1254  SCIP_ROUNDMODE roundmode;
1255  SCIP_INTERVAL prod;
1256  int i;
1257 
1258  roundmode = intervalGetRoundingMode();
1259 
1260  resultant->inf = 0.0;
1261  resultant->sup = 0.0;
1262 
1263  /* compute infimum of resultant */
1265  for( i = 0; i < length && resultant->inf > -infinity; ++i )
1266  {
1267  SCIPintervalSetEntire(infinity, &prod);
1268  SCIPintervalMulInf(infinity, &prod, operand1[i], operand2[i]);
1269  SCIPintervalAddInf(infinity, resultant, *resultant, prod);
1270  }
1271  assert(resultant->sup == 0.0);
1272 
1273  /* compute supremum of resultant */
1275  for( i = 0; i < length && resultant->sup < infinity ; ++i )
1276  {
1277  SCIPintervalSetEntire(infinity, &prod);
1278  SCIPintervalMulSup(infinity, &prod, operand1[i], operand2[i]);
1279  SCIPintervalAddSup(infinity, resultant, *resultant, prod);
1280  }
1281 
1282  intervalSetRoundingMode(roundmode);
1283 }
1284 
1285 /** computes scalar product of a vector of intervals and a vector of scalars and stores infimum of result in infimum of
1286  * resultant
1287  */
1289  SCIP_Real infinity, /**< value for infinity */
1290  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1291  int length, /**< length of vectors */
1292  SCIP_INTERVAL* operand1, /**< first vector as array of intervals */
1293  SCIP_Real* operand2 /**< second vector as array of scalars; can have +/-inf entries */
1294  )
1295 {
1296  SCIP_INTERVAL prod;
1297  int i;
1298 
1300 
1301  resultant->inf = 0.0;
1302 
1303  /* compute infimum of resultant */
1304  SCIPintervalSetEntire(infinity, &prod);
1305  for( i = 0; i < length && resultant->inf > -infinity; ++i )
1306  {
1307  SCIPintervalMulScalarInf(infinity, &prod, operand1[i], operand2[i]);
1308  assert(prod.sup >= infinity);
1309  SCIPintervalAddInf(infinity, resultant, *resultant, prod);
1310  }
1311 }
1312 
1313 /** computes the scalar product of a vector of intervals and a vector of scalars and stores supremum of result in
1314  * supremum of resultant
1315  */
1317  SCIP_Real infinity, /**< value for infinity */
1318  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1319  int length, /**< length of vectors */
1320  SCIP_INTERVAL* operand1, /**< first vector as array of intervals */
1321  SCIP_Real* operand2 /**< second vector as array of scalars; can have +/-inf entries */
1322  )
1323 {
1324  SCIP_INTERVAL prod;
1325  int i;
1326 
1328 
1329  resultant->sup = 0.0;
1330 
1331  /* compute supremum of resultant */
1332  SCIPintervalSetEntire(infinity, &prod);
1333  for( i = 0; i < length && resultant->sup < infinity; ++i )
1334  {
1335  SCIPintervalMulScalarSup(infinity, &prod, operand1[i], operand2[i]);
1336  assert(prod.inf <= -infinity);
1337  SCIPintervalAddSup(infinity, resultant, *resultant, prod);
1338  }
1339 }
1340 
1341 /** computes the scalar product of a vector of intervals and a vector of scalars and stores result in resultant */
1343  SCIP_Real infinity, /**< value for infinity */
1344  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1345  int length, /**< length of vectors */
1346  SCIP_INTERVAL* operand1, /**< first vector as array of intervals */
1347  SCIP_Real* operand2 /**< second vector as array of scalars; can have +/-inf entries */
1348  )
1349 {
1350  SCIP_ROUNDMODE roundmode;
1351 
1352  roundmode = intervalGetRoundingMode();
1353 
1354  resultant->inf = 0.0;
1355  resultant->sup = 0.0;
1356 
1357  /* compute infimum of resultant */
1359  SCIPintervalScalprodScalarsInf(infinity, resultant, length, operand1, operand2);
1360  assert(resultant->sup == 0.0);
1361 
1362  /* compute supremum of resultant */
1364  SCIPintervalScalprodScalarsSup(infinity, resultant, length, operand1, operand2);
1365 
1366  intervalSetRoundingMode(roundmode);
1367 }
1368 
1369 /** squares operand and stores result in resultant */
1371  SCIP_Real infinity, /**< value for infinity */
1372  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1373  SCIP_INTERVAL operand /**< operand of operation */
1374  )
1375 {
1376  SCIP_ROUNDMODE roundmode;
1377 
1378  assert(resultant != NULL);
1379  assert(!SCIPintervalIsEmpty(infinity, operand));
1380 
1381  roundmode = intervalGetRoundingMode();
1382 
1383  if( operand.sup <= 0.0 )
1384  { /* operand is left of 0.0 */
1385  if( operand.sup <= -infinity )
1386  resultant->inf = infinity;
1387  else
1388  {
1390  resultant->inf = operand.sup * operand.sup;
1391  }
1392 
1393  if( operand.inf <= -infinity )
1394  resultant->sup = infinity;
1395  else
1396  {
1398  resultant->sup = operand.inf * operand.inf;
1399  }
1400  }
1401  else if( operand.inf >= 0.0 )
1402  { /* operand is right of 0.0 */
1403  if( operand.inf >= infinity )
1404  resultant->inf = infinity;
1405  else
1406  {
1408  resultant->inf = operand.inf * operand.inf;
1409  }
1410 
1411  if( operand.sup >= infinity )
1412  resultant->sup = infinity;
1413  else
1414  {
1416  resultant->sup = operand.sup * operand.sup;
1417  }
1418  }
1419  else
1420  { /* [-,+]^2 */
1421  resultant->inf = 0.0;
1422  if( operand.inf <= -infinity || operand.sup >= infinity )
1423  resultant->sup = infinity;
1424  else
1425  {
1426  SCIP_Real x;
1427  SCIP_Real y;
1428 
1430  x = operand.inf * operand.inf;
1431  y = operand.sup * operand.sup;
1432  resultant->sup = MAX(x, y);
1433  }
1434  }
1435 
1436  intervalSetRoundingMode(roundmode);
1437 }
1438 
1439 /** stores (positive part of) square root of operand in resultant
1440  * @attention we assume a correctly rounded sqrt(double) function when rounding is to nearest
1441  */
1443  SCIP_Real infinity, /**< value for infinity */
1444  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1445  SCIP_INTERVAL operand /**< operand of operation */
1446  )
1447 {
1448  assert(resultant != NULL);
1449  assert(!SCIPintervalIsEmpty(infinity, operand));
1450 
1451  if( operand.sup < 0.0 )
1452  {
1453  SCIPintervalSetEmpty(resultant);
1454  return;
1455  }
1456 
1457  if( operand.inf == operand.sup ) /*lint !e777 */
1458  {
1459  if( operand.inf >= infinity )
1460  {
1461  resultant->inf = infinity;
1462  resultant->sup = infinity;
1463  }
1464  else
1465  {
1466  SCIP_Real tmp;
1467 
1468  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1469  tmp = sqrt(operand.inf);
1470  resultant->inf = SCIPnextafter(tmp, SCIP_REAL_MIN);
1471  resultant->sup = SCIPnextafter(tmp, SCIP_REAL_MAX);
1472  }
1473 
1474  return;
1475  }
1476 
1477  if( operand.inf <= 0.0 )
1478  resultant->inf = 0.0;
1479  else if( operand.inf >= infinity )
1480  {
1481  resultant->inf = infinity;
1482  resultant->sup = infinity;
1483  }
1484  else
1485  {
1486  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1487  resultant->inf = SCIPnextafter(sqrt(operand.inf), SCIP_REAL_MIN);
1488  }
1489 
1490  if( operand.sup >= infinity )
1491  resultant->sup = infinity;
1492  else
1493  {
1494  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1495  resultant->sup = SCIPnextafter(sqrt(operand.sup), SCIP_REAL_MAX);
1496  }
1497 }
1498 
1499 /** stores operand1 to the power of operand2 in resultant
1500  *
1501  * uses SCIPintervalPowerScalar if operand2 is a scalar, otherwise computes exp(op2*log(op1))
1502  */
1504  SCIP_Real infinity, /**< value for infinity */
1505  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1506  SCIP_INTERVAL operand1, /**< first operand of operation */
1507  SCIP_INTERVAL operand2 /**< second operand of operation */
1508  )
1509 { /*lint --e{777}*/
1510  assert(resultant != NULL);
1511  assert(!SCIPintervalIsEmpty(infinity, operand1));
1512  assert(!SCIPintervalIsEmpty(infinity, operand2));
1513 
1514  if( operand2.inf == operand2.sup )
1515  { /* operand is number */
1516  SCIPintervalPowerScalar(infinity, resultant, operand1, operand2.inf);
1517  return;
1518  }
1519 
1520  /* log([..,0]) will give an empty interval below, but we want [0,0]^exponent to be 0
1521  * if 0 is in exponent, then resultant should also contain 1 (the case exponent == [0,0] is handled above)
1522  */
1523  if( operand1.sup == 0.0 )
1524  {
1525  if( operand2.inf <= 0.0 && operand2.sup >= 0.0 )
1526  SCIPintervalSetBounds(resultant, 0.0, 1.0);
1527  else
1528  SCIPintervalSet(resultant, 0.0);
1529  return;
1530  }
1531 
1532  /* resultant := log(op1) */
1533  SCIPintervalLog(infinity, resultant, operand1);
1534  if( SCIPintervalIsEmpty(infinity, *resultant) )
1535  return;
1536 
1537  /* resultant := op2 * resultant */
1538  SCIPintervalMul(infinity, resultant, operand2, *resultant);
1539 
1540  /* resultant := exp(resultant) */
1541  SCIPintervalExp(infinity, resultant, *resultant);
1542 }
1543 
1544 /** computes lower bound on power of a scalar operand1 to an integer operand2
1545  * both operands need to be finite numbers
1546  * need to have operand1 >= 0 and need to have operand2 >= 0 if operand1 == 0
1547  */
1549  SCIP_Real operand1, /**< first operand of operation */
1550  int operand2 /**< second operand of operation */
1551  )
1552 {
1553  SCIP_ROUNDMODE roundmode;
1554  SCIP_Real result;
1555 
1556  assert(operand1 >= 0.0);
1557 
1558  if( operand1 == 0.0 )
1559  {
1560  assert(operand2 >= 0);
1561  if( operand2 == 0 )
1562  return 1.0; /* 0^0 = 1 */
1563  else
1564  return 0.0; /* 0^positive = 0 */
1565  }
1566 
1567  /* 1^n = 1, x^0 = 1 */
1568  if( operand1 == 1.0 || operand2 == 0 )
1569  return 1.0;
1570 
1571  if( operand2 < 0 )
1572  {
1573  /* x^n = 1 / x^(-n) */
1574  result = SCIPintervalPowerScalarIntegerSup(operand1, -operand2);
1575  assert(result != 0.0);
1576 
1577  roundmode = intervalGetRoundingMode();
1579  result = 1.0 / result;
1580  intervalSetRoundingMode(roundmode);
1581  }
1582  else
1583  {
1584  unsigned int n;
1585  SCIP_Real z;
1586 
1587  roundmode = intervalGetRoundingMode();
1588 
1589  result = 1.0;
1590  n = (unsigned int)operand2;
1591  z = operand1;
1592 
1594 
1595  /* use a binary exponentiation algorithm:
1596  * consider the binary representation of n: n = sum_i 2^i d_i with d_i \in {0,1}
1597  * then x^n = prod_{i:d_i=1} x^(2^i)
1598  * in the following, we loop over i=1,..., thereby storing x^(2^i) in z
1599  * whenever d_i is 1, we multiply result with x^(2^i) (the current z)
1600  * at the last (highest) i with d_i = 1 we stop, thus having x^n stored in result
1601  *
1602  * the binary representation of n and bit shifting is used for the loop
1603  */
1604  assert(n >= 1);
1605  do
1606  {
1607  if( n & 1 ) /* n is odd (d_i=1), so multiply result with current z (=x^{2^i}) */
1608  {
1609  result = result * z;
1610  n >>= 1;
1611  if( n == 0 )
1612  break;
1613  }
1614  else
1615  n >>= 1;
1616  z = z * z;
1617  }
1618  while( TRUE ); /*lint !e506 */
1619 
1620  intervalSetRoundingMode(roundmode);
1621  }
1622 
1623  return result;
1624 }
1625 
1626 /** computes upper bound on power of a scalar operand1 to an integer operand2
1627  * both operands need to be finite numbers
1628  * need to have operand1 >= 0 and need to have operand2 >= 0 if operand1 == 0
1629  */
1631  SCIP_Real operand1, /**< first operand of operation */
1632  int operand2 /**< second operand of operation */
1633  )
1634 {
1635  SCIP_ROUNDMODE roundmode;
1636  SCIP_Real result;
1637 
1638  assert(operand1 >= 0.0);
1639 
1640  if( operand1 == 0.0 )
1641  {
1642  assert(operand2 >= 0);
1643  if( operand2 == 0 )
1644  return 1.0; /* 0^0 = 1 */
1645  else
1646  return 0.0; /* 0^positive = 0 */
1647  }
1648 
1649  /* 1^n = 1, x^0 = 1 */
1650  if( operand1 == 1.0 || operand2 == 0 )
1651  return 1.0;
1652 
1653  if( operand2 < 0 )
1654  {
1655  /* x^n = 1 / x^(-n) */
1656  result = SCIPintervalPowerScalarIntegerInf(operand1, -operand2);
1657  assert(result != 0.0);
1658 
1659  roundmode = intervalGetRoundingMode();
1661  result = 1.0 / result;
1662  intervalSetRoundingMode(roundmode);
1663  }
1664  else
1665  {
1666  unsigned int n;
1667  SCIP_Real z;
1668 
1669  roundmode = intervalGetRoundingMode();
1670 
1671  result = 1.0;
1672  n = (unsigned int)operand2;
1673  z = operand1;
1674 
1676 
1677  /* use a binary exponentiation algorithm... see comments in SCIPintervalPowerScalarIntegerInf */
1678  assert(n >= 1);
1679  do
1680  {
1681  if( n&1 )
1682  {
1683  result = result * z;
1684  n >>= 1;
1685  if( n == 0 )
1686  break;
1687  }
1688  else
1689  n >>= 1;
1690  z = z * z;
1691  }
1692  while( TRUE ); /*lint !e506 */
1693 
1694  intervalSetRoundingMode(roundmode);
1695  }
1696 
1697  return result;
1698 }
1699 
1700 /** computes bounds on power of a scalar operand1 to an integer operand2
1701  * both operands need to be finite numbers
1702  * need to have operand1 >= 0 and need to have operand2 >= 0 if operand1 == 0
1703  */
1705  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1706  SCIP_Real operand1, /**< first operand of operation */
1707  int operand2 /**< second operand of operation */
1708  )
1709 {
1710  SCIP_ROUNDMODE roundmode;
1711 
1712  assert(operand1 >= 0.0);
1713 
1714  if( operand1 == 0.0 )
1715  {
1716  assert(operand2 >= 0);
1717  if( operand2 == 0 )
1718  {
1719  SCIPintervalSet(resultant, 1.0); /* 0^0 = 1 */
1720  return;
1721  }
1722  else
1723  {
1724  SCIPintervalSet(resultant, 0.0); /* 0^positive = 0 */
1725  return;
1726  }
1727  }
1728 
1729  /* 1^n = 1, x^0 = 1 */
1730  if( operand1 == 1.0 || operand2 == 0 )
1731  {
1732  SCIPintervalSet(resultant, 1.0);
1733  return;
1734  }
1735 
1736  if( operand2 < 0 )
1737  {
1738  /* x^n = 1 / x^(-n) */
1739  SCIPintervalPowerScalarInteger(resultant, operand1, -operand2);
1740  assert(resultant->inf > 0.0 || resultant->sup < 0.0);
1741  SCIPintervalReciprocal(SCIP_REAL_MAX, resultant, *resultant); /* value for infinity does not matter, since there should be no 0.0 in the interval, so just use something large enough */
1742  }
1743  else
1744  {
1745  unsigned int n;
1746  SCIP_Real z_inf;
1747  SCIP_Real z_sup;
1748  SCIP_Real result_sup;
1749  SCIP_Real result_inf;
1750 
1751  roundmode = intervalGetRoundingMode();
1752 
1753  result_inf = 1.0;
1754  result_sup = 1.0;
1755  z_inf = operand1;
1756  z_sup = operand1;
1757  n = (unsigned int)operand2;
1758 
1760 
1761  /* use a binary exponentiation algorithm... see comments in SCIPintervalPowerScalarIntegerInf
1762  * we compute lower and upper bounds within the same loop
1763  * to get correct lower bounds while rounding mode is upwards, we negate arguments */
1764  assert(n >= 1);
1765  do
1766  {
1767  if( n & 1 )
1768  {
1769  result_inf = negate(negate(result_inf) * z_inf);
1770  result_sup = result_sup * z_sup;
1771  n >>= 1;
1772  if( n == 0 )
1773  break;
1774  }
1775  else
1776  n >>= 1;
1777  z_inf = negate(negate(z_inf) * z_inf);
1778  z_sup = z_sup * z_sup;
1779  }
1780  while( TRUE ); /*lint !e506 */
1781 
1782  intervalSetRoundingMode(roundmode);
1783 
1784  resultant->inf = result_inf;
1785  resultant->sup = result_sup;
1786  assert(resultant->inf <= resultant->sup);
1787  }
1788 }
1789 
1790 /** stores bounds on the power of a scalar operand1 to a scalar operand2 in resultant
1791  * both operands need to be finite numbers
1792  * need to have operand1 >= 0 or operand2 integer and need to have operand2 >= 0 if operand1 == 0
1793  * @attention we assume a correctly rounded pow(double) function when rounding is to nearest
1794  */
1796  SCIP_INTERVAL* resultant, /**< resultant of operation */
1797  SCIP_Real operand1, /**< first operand of operation */
1798  SCIP_Real operand2 /**< second operand of operation */
1799  )
1800 {
1801  SCIP_Real result;
1802 
1803  assert(resultant != NULL);
1804 
1805  if( operand1 == 0.0 )
1806  {
1807  assert(operand2 >= 0);
1808  if( operand2 == 0 )
1809  {
1810  SCIPintervalSet(resultant, 1.0); /* 0^0 = 1 */
1811  return;
1812  }
1813  else
1814  {
1815  SCIPintervalSet(resultant, 0.0); /* 0^positive = 0 */
1816  return;
1817  }
1818  }
1819 
1820  /* 1^n = 1, x^0 = 1 */
1821  if( operand1 == 1.0 || operand2 == 0 )
1822  {
1823  SCIPintervalSet(resultant, 1.0);
1824  return;
1825  }
1826 
1827  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1828  result = pow(operand1, operand2);
1829 
1830  /* to get safe bounds, get the floating point numbers just below and above result */
1831  resultant->inf = SCIPnextafter(result, SCIP_REAL_MIN);
1832  resultant->sup = SCIPnextafter(result, SCIP_REAL_MAX);
1833 }
1834 
1835 /** stores operand1 to the power of the scalar operand2 in resultant
1836  * @attention we assume a correctly rounded pow(double) function when rounding is to nearest
1837  */
1839  SCIP_Real infinity, /**< value for infinity */
1840  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1841  SCIP_INTERVAL operand1, /**< first operand of operation */
1842  SCIP_Real operand2 /**< second operand of operation */
1843  )
1844 { /*lint --e{777}*/
1845  SCIP_Bool op2isint;
1846 
1847  assert(resultant != NULL);
1848  assert(!SCIPintervalIsEmpty(infinity, operand1));
1849 
1850  if( operand2 == infinity )
1851  {
1852  /* 0^infinity = 0
1853  * +^infinity = infinity
1854  * -^infinity = -infinity
1855  */
1856  if( operand1.inf < 0.0 )
1857  resultant->inf = -infinity;
1858  else
1859  resultant->inf = 0.0;
1860  if( operand1.sup > 0.0 )
1861  resultant->sup = infinity;
1862  else
1863  resultant->sup = 0.0;
1864  return;
1865  }
1866 
1867  if( operand2 == 0.0 )
1868  { /* special case, since x^0 = 1 for x != 0, but 0^0 = 0 */
1869  if( operand1.inf == 0.0 && operand1.sup == 0.0 )
1870  {
1871  resultant->inf = 0.0;
1872  resultant->sup = 0.0;
1873  }
1874  else if( operand1.inf <= 0.0 || operand1.sup >= 0.0 )
1875  { /* 0.0 in x gives [0,1] */
1876  resultant->inf = 0.0;
1877  resultant->sup = 1.0;
1878  }
1879  else
1880  { /* 0.0 outside x gives [1,1] */
1881  resultant->inf = 1.0;
1882  resultant->sup = 1.0;
1883  }
1884  return;
1885  }
1886 
1887  if( operand2 == 1.0 )
1888  {
1889  /* x^1 = x */
1890  *resultant = operand1;
1891  return;
1892  }
1893 
1894  op2isint = (ceil(operand2) == operand2);
1895 
1896  if( !op2isint && operand1.inf < 0.0 )
1897  { /* x^n with x negative not defined for n not integer*/
1898  operand1.inf = 0.0;
1899  if( operand1.sup < operand1.inf )
1900  {
1901  SCIPintervalSetEmpty(resultant);
1902  return;
1903  }
1904  }
1905 
1906  if( operand1.inf >= 0.0 )
1907  { /* easy case: x^n with x >= 0 */
1908  if( operand2 >= 0.0 )
1909  {
1910  /* inf^+ = inf */
1911  if( operand1.inf >= infinity )
1912  resultant->inf = infinity;
1913  else if( operand1.inf > 0.0 )
1914  {
1915  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1916  resultant->inf = SCIPnextafter(pow(operand1.inf, operand2), SCIP_REAL_MIN);
1917  }
1918  else
1919  resultant->inf = 0.0;
1920 
1921  if( operand1.sup >= infinity )
1922  resultant->sup = infinity;
1923  else if( operand1.sup > 0.0 )
1924  {
1925  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1926  resultant->sup = SCIPnextafter(pow(operand1.sup, operand2), SCIP_REAL_MAX);
1927  }
1928  else
1929  resultant->sup = 0.0;
1930  }
1931  else
1932  {
1933  if( operand1.sup >= infinity )
1934  resultant->inf = 0.0;
1935  else if( operand1.sup == 0.0 )
1936  {
1937  /* x^(negative even) = infinity for x->0 (from both sides),
1938  * but x^(negative odd) = -infinity for x->0 from left side */
1939  if( ceil(operand2/2) == operand2/2 )
1940  resultant->inf = infinity;
1941  else
1942  resultant->inf = -infinity;
1943  }
1944  else
1945  {
1946  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1947  resultant->inf = SCIPnextafter(pow(operand1.sup, operand2), SCIP_REAL_MIN);
1948  }
1949 
1950  /* 0^(negative) = infinity */
1951  if( operand1.inf == 0.0 )
1952  resultant->sup = infinity;
1953  else
1954  {
1955  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1956  resultant->sup = SCIPnextafter(pow(operand1.inf, operand2), SCIP_REAL_MAX);
1957  }
1958  }
1959  }
1960  else if( operand1.sup <= 0.0 )
1961  { /* more difficult case: x^n with x < 0; we now know, that n is integer */
1962  assert(op2isint);
1963  if( operand2 >= 0.0 && ceil(operand2/2) == operand2/2 )
1964  {
1965  /* x^n with n>=2 and even -> x^n is monotonically decreasing for x < 0 */
1966  if( operand1.sup == -infinity )
1967  /* (-inf)^n = inf */
1968  resultant->inf = infinity;
1969  else
1970  resultant->inf = SCIPintervalPowerScalarIntegerInf(-operand1.sup, (int)operand2);
1971 
1972  if( operand1.inf <= -infinity )
1973  resultant->sup = infinity;
1974  else
1975  resultant->sup = SCIPintervalPowerScalarIntegerSup(-operand1.inf, (int)operand2);
1976  }
1977  else if( operand2 <= 0.0 && ceil(operand2/2) != operand2/2 )
1978  {
1979  /* x^n with n<=-1 and odd -> x^n = 1/x^(-n) is monotonically decreasing for x<0 */
1980  if( operand1.sup == -infinity )
1981  /* (-inf)^n = 1/(-inf)^(-n) = 1/(-inf) = 0 */
1982  resultant->inf = 0.0;
1983  else if( operand1.sup == 0.0 )
1984  /* x^n -> -infinity for x->0 from left */
1985  resultant->inf = -infinity;
1986  else
1987  resultant->inf = -SCIPintervalPowerScalarIntegerSup(-operand1.sup, (int)operand2);
1988 
1989  if( operand1.inf <= -infinity )
1990  /* (-inf)^n = 1/(-inf)^(-n) = 1/(-inf) = 0 */
1991  resultant->sup = 0.0;
1992  else if( operand1.inf == 0.0 )
1993  /* x^n -> infinity for x->0 from right */
1994  resultant->sup = infinity;
1995  else
1996  resultant->sup = -SCIPintervalPowerScalarIntegerInf(-operand1.inf, (int)operand2);
1997  }
1998  else if( operand2 >= 0.0 )
1999  {
2000  /* x^n with n>0 and odd -> x^n is monotonically increasing for x<0 */
2001  if( operand1.inf <= -infinity )
2002  resultant->inf = -infinity;
2003  else
2004  resultant->inf = -SCIPintervalPowerScalarIntegerSup(-operand1.inf, (int)operand2);
2005 
2006  if( operand1.sup <= -infinity )
2007  resultant->sup = -infinity;
2008  else
2009  resultant->sup = -SCIPintervalPowerScalarIntegerInf(-operand1.sup, (int)operand2);
2010  }
2011  else
2012  {
2013  /* x^n with n<0 and even -> x^n is monotonically increasing for x<0 */
2014  if( operand1.inf <= -infinity )
2015  resultant->inf = 0.0;
2016  else if( operand1.inf == 0.0 )
2017  /* x^n -> infinity for x->0 from both sides */
2018  resultant->inf = infinity;
2019  else
2020  resultant->inf = SCIPintervalPowerScalarIntegerSup(-operand1.inf, (int)operand2);
2021 
2022  if( operand1.sup <= -infinity )
2023  resultant->sup = 0.0;
2024  else if( operand1.sup == 0.0 )
2025  /* x^n -> infinity for x->0 from both sides */
2026  resultant->sup = infinity;
2027  else
2028  resultant->sup = SCIPintervalPowerScalarIntegerSup(-operand1.sup, (int)operand2);
2029  }
2030  assert(resultant->inf <= resultant->sup || resultant->inf >= infinity || resultant->sup <= -infinity);
2031  }
2032  else
2033  { /* similar difficult case: x^n with x in [<0, >0], but n is integer */
2034  assert(op2isint); /* otherwise we had set operand1.inf == 0.0, which was handled in first case */
2035  if( operand2 >= 0.0 && operand2/2 == ceil(operand2/2) )
2036  {
2037  /* n even positive integer */
2038  resultant->inf = 0.0;
2039  if( operand1.inf == -infinity || operand1.sup == infinity )
2040  resultant->sup = infinity;
2041  else
2042  resultant->sup = SCIPintervalPowerScalarIntegerSup(MAX(-operand1.inf, operand1.sup), (int)operand2);
2043  }
2044  else if( operand2 <= 0.0 && ceil(operand2/2) == operand2/2 )
2045  {
2046  /* n even negative integer */
2047  resultant->sup = infinity; /* since 0^n = infinity */
2048  if( operand1.inf == -infinity || operand1.sup == infinity )
2049  resultant->inf = 0.0;
2050  else
2051  resultant->inf = SCIPintervalPowerScalarIntegerInf(MAX(-operand1.inf, operand1.sup), (int)operand2);
2052  }
2053  else if( operand2 >= 0.0 )
2054  {
2055  /* n odd positive integer, so monotonically increasing function */
2056  if( operand1.inf == -infinity )
2057  resultant->inf = -infinity;
2058  else
2059  resultant->inf = -SCIPintervalPowerScalarIntegerSup(-operand1.inf, (int)operand2);
2060  if( operand1.sup == infinity )
2061  resultant->sup = infinity;
2062  else
2063  resultant->sup = SCIPintervalPowerScalarIntegerSup(operand1.sup, (int)operand2);
2064  }
2065  else
2066  {
2067  /* n odd negative integer:
2068  * x^n -> -infinity for x->0 from left
2069  * x^n -> infinity for x->0 from right */
2070  resultant->inf = -infinity;
2071  resultant->sup = infinity;
2072  }
2073  }
2074 
2075  /* if value for infinity is too small, relax intervals so they do not appear empty */
2076  if( resultant->inf > infinity )
2077  resultant->inf = infinity;
2078  if( resultant->sup < -infinity )
2079  resultant->sup = -infinity;
2080 }
2081 
2082 /** given an interval for the image of a power operation, computes an interval for the origin
2083  * that is, for y = x^p with p = exponent a given scalar and y = image a given interval,
2084  * computes a subinterval x of basedomain such that y in x^p and such that for all z in basedomain less x, z^p not in y
2085  */
2087  SCIP_Real infinity, /**< value for infinity */
2088  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2089  SCIP_INTERVAL basedomain, /**< domain of base */
2090  SCIP_Real exponent, /**< exponent */
2091  SCIP_INTERVAL image /**< interval image of power */
2092  )
2093 {
2094  SCIP_INTERVAL tmp;
2095  SCIP_INTERVAL exprecip;
2096 
2097  assert(resultant != NULL);
2098  assert(image.inf <= image.sup);
2099  assert(basedomain.inf <= basedomain.sup);
2100 
2101  if( exponent == 0.0 )
2102  {
2103  /* exponent is 0.0 */
2104  if( image.inf <= 1.0 && image.sup >= 1.0 )
2105  {
2106  /* 1 in image -> resultant = entire */
2107  *resultant = basedomain;
2108  }
2109  else if( image.inf <= 0.0 && image.sup >= 0.0 )
2110  {
2111  /* 0 in image, 1 not in image -> resultant = 0 (provided 0^0 = 0 ???)
2112  * -> resultant = {0} intersected with basedomain */
2113  SCIPintervalSetBounds(resultant, MAX(0.0, basedomain.inf), MIN(0.0, basedomain.sup));
2114  }
2115  else
2116  {
2117  /* 0 and 1 not in image -> resultant = empty */
2118  SCIPintervalSetEmpty(resultant);
2119  }
2120  return;
2121  }
2122 
2123  /* i = b^e
2124  * i >= 0 -> b = i^(1/e) [union -i^(1/e), if e is even]
2125  * i < 0, e odd integer -> b = -(-i)^(1/e)
2126  * i < 0, e even integer or fractional -> empty
2127  */
2128 
2129  SCIPintervalSetBounds(&exprecip, exponent, exponent);
2130  SCIPintervalReciprocal(infinity, &exprecip, exprecip);
2131 
2132  /* invert positive part of image, if any */
2133  if( image.sup >= 0.0 )
2134  {
2135  SCIPintervalSetBounds(&tmp, MAX(image.inf, 0.0), image.sup);
2136  SCIPintervalPower(infinity, resultant, tmp, exprecip);
2137  if( basedomain.inf <= -resultant->inf && EPSISINT(exponent, 0.0) && (int)exponent % 2 == 0 ) /*lint !e835 */
2138  {
2139  if( basedomain.sup < resultant->inf )
2140  SCIPintervalSetBounds(resultant, -resultant->sup, -resultant->inf);
2141  else
2142  SCIPintervalSetBounds(resultant, -resultant->sup, resultant->sup);
2143  }
2144 
2145  SCIPintervalIntersect(resultant, *resultant, basedomain);
2146  }
2147  else
2148  SCIPintervalSetEmpty(resultant);
2149 
2150  /* invert negative part of image, if any and if base can take negative value and if exponent is such that negative values are possible */
2151  if( image.inf < 0.0 && basedomain.inf < 0.0 && EPSISINT(exponent, 0.0) && ((int)exponent % 2 != 0) ) /*lint !e835 */
2152  {
2153  SCIPintervalSetBounds(&tmp, MAX(-image.sup, 0.0), -image.inf);
2154  SCIPintervalPower(infinity, &tmp, tmp, exprecip);
2155  SCIPintervalSetBounds(&tmp, -tmp.sup, -tmp.inf);
2156  SCIPintervalIntersect(&tmp, basedomain, tmp);
2157  SCIPintervalUnify(resultant, *resultant, tmp);
2158  }
2159 }
2160 
2161 /** stores operand1 to the signed power of the scalar positive operand2 in resultant
2162  *
2163  * the signed power of x w.r.t. an exponent n >= 0 is given as sign(x) * abs(x)^n
2164  *
2165  * @attention we assume correctly rounded sqrt(double) and pow(double) functions when rounding is to nearest
2166  */
2168  SCIP_Real infinity, /**< value for infinity */
2169  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2170  SCIP_INTERVAL operand1, /**< first operand of operation */
2171  SCIP_Real operand2 /**< second operand of operation */
2172  )
2173 {
2174  SCIP_ROUNDMODE roundmode;
2175  assert(resultant != NULL);
2176 
2177  assert(!SCIPintervalIsEmpty(infinity, operand1));
2178  assert(operand2 >= 0.0);
2179 
2180  if( operand2 == infinity ) /*lint !e777 */
2181  {
2182  /* 0^infinity = 0
2183  * +^infinity = infinity
2184  *-+^infinity = -infinity
2185  */
2186  if( operand1.inf < 0.0 )
2187  resultant->inf = -infinity;
2188  else
2189  resultant->inf = 0.0;
2190  if( operand1.sup > 0.0 )
2191  resultant->sup = infinity;
2192  else
2193  resultant->sup = 0.0;
2194  return;
2195  }
2196 
2197  if( operand2 == 0.0 )
2198  {
2199  /* special case, since x^0 = 1 for x != 0, but 0^0 = 0 */
2200  if( operand1.inf < 0.0 )
2201  resultant->inf = -1.0;
2202  else if( operand1.inf == 0.0 )
2203  resultant->inf = 0.0;
2204  else
2205  resultant->inf = 1.0;
2206 
2207  if( operand1.sup < 0.0 )
2208  resultant->sup = -1.0;
2209  else if( operand1.sup == 0.0 )
2210  resultant->sup = 0.0;
2211  else
2212  resultant->sup = 1.0;
2213 
2214  return;
2215  }
2216 
2217  if( operand2 == 1.0 )
2218  { /* easy case that should run fast */
2219  *resultant = operand1;
2220  return;
2221  }
2222 
2223  roundmode = intervalGetRoundingMode();
2224 
2225  if( operand2 == 2.0 )
2226  { /* common case where pow can easily be avoided */
2227  if( operand1.inf <= -infinity )
2228  {
2229  resultant->inf = -infinity;
2230  }
2231  else if( operand1.inf >= infinity )
2232  {
2233  resultant->inf = infinity;
2234  }
2235  else if( operand1.inf > 0.0 )
2236  {
2238  resultant->inf = operand1.inf * operand1.inf;
2239  }
2240  else
2241  {
2242  /* need upwards since we negate result of multiplication */
2244  resultant->inf = negate(operand1.inf * operand1.inf);
2245  }
2246 
2247  if( operand1.sup >= infinity )
2248  {
2249  resultant->sup = infinity;
2250  }
2251  else if( operand1.sup <= -infinity )
2252  {
2253  resultant->sup = -infinity;
2254  }
2255  else if( operand1.sup > 0.0 )
2256  {
2258  resultant->sup = operand1.sup * operand1.sup;
2259  }
2260  else
2261  {
2262  /* need downwards since we negate result of multiplication */
2264  resultant->sup = negate(operand1.sup * operand1.sup);
2265  }
2266  assert(resultant->inf <= resultant->sup);
2267  }
2268  else if( operand2 == 0.5 )
2269  { /* another common case where pow can easily be avoided */
2270  if( operand1.inf <= -infinity )
2271  resultant->inf = -infinity;
2272  else if( operand1.inf >= infinity )
2273  resultant->inf = infinity;
2274  else if( operand1.inf >= 0.0 )
2275  {
2277  resultant->inf = SCIPnextafter(sqrt( operand1.inf), SCIP_REAL_MIN);
2278  }
2279  else
2280  {
2282  resultant->inf = -SCIPnextafter(sqrt(-operand1.inf), SCIP_REAL_MAX);
2283  }
2284 
2285  if( operand1.sup >= infinity )
2286  resultant->sup = infinity;
2287  else if( operand1.sup <= -infinity )
2288  resultant->sup = -infinity;
2289  else if( operand1.sup > 0.0 )
2290  {
2292  resultant->sup = SCIPnextafter(sqrt( operand1.sup), SCIP_REAL_MAX);
2293  }
2294  else
2295  {
2297  resultant->sup = -SCIPnextafter(sqrt(-operand1.sup), SCIP_REAL_MAX);
2298  }
2299  assert(resultant->inf <= resultant->sup);
2300  }
2301  else
2302  {
2303  if( operand1.inf <= -infinity )
2304  resultant->inf = -infinity;
2305  else if( operand1.inf >= infinity )
2306  resultant->inf = infinity;
2307  else if( operand1.inf > 0.0 )
2308  {
2310  resultant->inf = SCIPnextafter(pow( operand1.inf, operand2), SCIP_REAL_MIN);
2311  }
2312  else
2313  {
2315  resultant->inf = -SCIPnextafter(pow(-operand1.inf, operand2), SCIP_REAL_MAX);
2316  }
2317 
2318  if( operand1.sup >= infinity )
2319  resultant->sup = infinity;
2320  else if( operand1.sup <= -infinity )
2321  resultant->sup = -infinity;
2322  else if( operand1.sup > 0.0 )
2323  {
2325  resultant->sup = SCIPnextafter(pow( operand1.sup, operand2), SCIP_REAL_MAX);
2326  }
2327  else
2328  {
2330  resultant->sup = -SCIPnextafter(pow(-operand1.sup, operand2), SCIP_REAL_MIN);
2331  }
2332  }
2333 
2334  intervalSetRoundingMode(roundmode);
2335 }
2336 
2337 /** computes the reciprocal of an interval
2338  */
2340  SCIP_Real infinity, /**< value for infinity */
2341  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2342  SCIP_INTERVAL operand /**< operand of operation */
2343  )
2344 {
2345  SCIP_ROUNDMODE roundmode;
2346 
2347  assert(resultant != NULL);
2348  assert(!SCIPintervalIsEmpty(infinity, operand));
2349 
2350  if( operand.inf == 0.0 && operand.sup == 0.0 )
2351  { /* 1/0 = [-inf,inf] */
2352  resultant->inf = infinity;
2353  resultant->sup = -infinity;
2354  return;
2355  }
2356 
2357  roundmode = intervalGetRoundingMode();
2358 
2359  if( operand.inf >= 0.0 )
2360  { /* 1/x with x >= 0 */
2361  if( operand.sup >= infinity )
2362  resultant->inf = 0.0;
2363  else
2364  {
2366  resultant->inf = 1.0 / operand.sup;
2367  }
2368 
2369  if( operand.inf >= infinity )
2370  resultant->sup = 0.0;
2371  else if( operand.inf == 0.0 )
2372  resultant->sup = infinity;
2373  else
2374  {
2376  resultant->sup = 1.0 / operand.inf;
2377  }
2378 
2379  intervalSetRoundingMode(roundmode);
2380  }
2381  else if( operand.sup <= 0.0 )
2382  { /* 1/x with x <= 0 */
2383  if( operand.sup <= -infinity )
2384  resultant->inf = 0.0;
2385  else if( operand.sup == 0.0 )
2386  resultant->inf = -infinity;
2387  else
2388  {
2390  resultant->inf = 1.0 / operand.sup;
2391  }
2392 
2393  if( operand.inf <= -infinity )
2394  resultant->sup = infinity;
2395  else
2396  {
2398  resultant->sup = 1.0 / operand.inf;
2399  }
2400  intervalSetRoundingMode(roundmode);
2401  }
2402  else
2403  { /* 1/x with x in [-,+] is division by zero */
2404  resultant->inf = -infinity;
2405  resultant->sup = infinity;
2406  }
2407 }
2408 
2409 /** stores exponential of operand in resultant
2410  * @attention we assume a correctly rounded exp(double) function when rounding is to nearest
2411  */
2413  SCIP_Real infinity, /**< value for infinity */
2414  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2415  SCIP_INTERVAL operand /**< operand of operation */
2416  )
2417 {
2418  assert(resultant != NULL);
2419  assert(!SCIPintervalIsEmpty(infinity, operand));
2420 
2421  if( operand.sup <= -infinity )
2422  {
2423  resultant->inf = 0.0;
2424  resultant->sup = 0.0;
2425  return;
2426  }
2427 
2428  if( operand.inf >= infinity )
2429  {
2430  resultant->inf = infinity;
2431  resultant->sup = infinity;
2432  return;
2433  }
2434 
2435  if( operand.inf == operand.sup ) /*lint !e777 */
2436  {
2437  if( operand.inf == 0.0 )
2438  {
2439  resultant->inf = 1.0;
2440  resultant->sup = 1.0;
2441  }
2442  else
2443  {
2444  SCIP_Real tmp;
2445 
2447  tmp = exp(operand.inf);
2448  resultant->inf = tmp > 0.0 ? SCIPnextafter(tmp, SCIP_REAL_MIN) : 0.0;
2449  assert(resultant->inf >= 0.0);
2450  resultant->sup = SCIPnextafter(tmp, SCIP_REAL_MAX);
2451 
2452  return;
2453  }
2454  }
2455 
2456  if( operand.inf <= -infinity )
2457  {
2458  resultant->inf = 0.0;
2459  }
2460  else if( operand.inf == 0.0 )
2461  {
2462  resultant->inf = 1.0;
2463  }
2464  else
2465  {
2466  SCIP_Real tmp;
2467 
2469  tmp = exp(operand.inf);
2470  resultant->inf = tmp > 0.0 ? SCIPnextafter(tmp, SCIP_REAL_MIN) : 0.0;
2471  /* make sure we do not exceed value for infinity, so interval is not declared as empty if inf and sup are both > infinity */
2472  if( resultant->inf >= infinity )
2473  resultant->inf = infinity;
2474  }
2475 
2476  if( operand.sup >= infinity )
2477  {
2478  resultant->sup = infinity;
2479  }
2480  else if( operand.sup == 0.0 )
2481  {
2482  resultant->sup = 1.0;
2483  }
2484  else
2485  {
2487  resultant->sup = SCIPnextafter(exp(operand.sup), SCIP_REAL_MAX);
2488  if( resultant->sup < -infinity )
2489  resultant->sup = -infinity;
2490  }
2491 }
2492 
2493 /** stores natural logarithm of operand in resultant
2494  * @attention we assume a correctly rounded log(double) function when rounding is to nearest
2495  */
2497  SCIP_Real infinity, /**< value for infinity */
2498  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2499  SCIP_INTERVAL operand /**< operand of operation */
2500  )
2501 {
2502  assert(resultant != NULL);
2503  assert(!SCIPintervalIsEmpty(infinity, operand));
2504 
2505  /* if operand.sup == 0.0, we could return -inf in resultant->sup, but that
2506  * seems of little use and just creates problems somewhere else, e.g., #1230
2507  */
2508  if( operand.sup <= 0.0 )
2509  {
2510  SCIPintervalSetEmpty(resultant);
2511  return;
2512  }
2513 
2514  if( operand.inf == operand.sup ) /*lint !e777 */
2515  {
2516  if( operand.sup == 1.0 )
2517  {
2518  resultant->inf = 0.0;
2519  resultant->sup = 0.0;
2520  }
2521  else
2522  {
2523  SCIP_Real tmp;
2524 
2526  tmp = log(operand.inf);
2527  resultant->inf = SCIPnextafter(tmp, SCIP_REAL_MIN);
2528  resultant->sup = SCIPnextafter(tmp, SCIP_REAL_MAX);
2529  }
2530 
2531  return;
2532  }
2533 
2534  if( operand.inf <= 0.0 )
2535  {
2536  resultant->inf = -infinity;
2537  }
2538  else if( operand.inf == 1.0 )
2539  {
2540  resultant->inf = 0.0;
2541  }
2542  else
2543  {
2545  resultant->inf = SCIPnextafter(log(operand.inf), SCIP_REAL_MIN);
2546  }
2547 
2548  if( operand.sup >= infinity )
2549  {
2550  resultant->sup = infinity;
2551  }
2552  else if( operand.sup == 1.0 )
2553  {
2554  resultant->sup = 0.0;
2555  }
2556  else
2557  {
2559  resultant->sup = SCIPnextafter(log(operand.sup), SCIP_REAL_MAX);
2560  }
2561 }
2562 
2563 /** stores minimum of operands in resultant */
2565  SCIP_Real infinity, /**< value for infinity */
2566  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2567  SCIP_INTERVAL operand1, /**< first operand of operation */
2568  SCIP_INTERVAL operand2 /**< second operand of operation */
2569  )
2570 {
2571  assert(resultant != NULL);
2572  assert(!SCIPintervalIsEmpty(infinity, operand1));
2573  assert(!SCIPintervalIsEmpty(infinity, operand2));
2574 
2575  resultant->inf = MIN(operand1.inf, operand2.inf);
2576  resultant->sup = MIN(operand1.sup, operand2.sup);
2577 }
2578 
2579 /** stores maximum of operands in resultant */
2581  SCIP_Real infinity, /**< value for infinity */
2582  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2583  SCIP_INTERVAL operand1, /**< first operand of operation */
2584  SCIP_INTERVAL operand2 /**< second operand of operation */
2585  )
2586 {
2587  assert(resultant != NULL);
2588  assert(!SCIPintervalIsEmpty(infinity, operand1));
2589  assert(!SCIPintervalIsEmpty(infinity, operand2));
2590 
2591  resultant->inf = MAX(operand1.inf, operand2.inf);
2592  resultant->sup = MAX(operand1.sup, operand2.sup);
2593 }
2594 
2595 /** stores absolute value of operand in resultant */
2597  SCIP_Real infinity, /**< value for infinity */
2598  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2599  SCIP_INTERVAL operand /**< operand of operation */
2600  )
2601 {
2602  assert(resultant != NULL);
2603  assert(!SCIPintervalIsEmpty(infinity, operand));
2604 
2605  if( operand.inf <= 0.0 && operand.sup >= 0.0)
2606  {
2607  resultant->inf = 0.0;
2608  resultant->sup = MAX(-operand.inf, operand.sup);
2609  }
2610  else if( operand.inf > 0.0 )
2611  {
2612  *resultant = operand;
2613  }
2614  else
2615  {
2616  resultant->inf = -operand.sup;
2617  resultant->sup = -operand.inf;
2618  }
2619 }
2620 
2621 /** stores sine value of operand in resultant
2622  * NOTE: the operations are not applied rounding-safe here
2623  */
2625  SCIP_Real infinity, /**< value for infinity */
2626  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2627  SCIP_INTERVAL operand /**< operand of operation */
2628  )
2629 {
2630  SCIP_Real intervallen;
2631  SCIP_Real modinf;
2632  SCIP_Real modsup;
2633  SCIP_Real finf;
2634  SCIP_Real fsup;
2635  int a;
2636  int b;
2637  int nbetween;
2638 
2639  assert(resultant != NULL);
2640  assert(!SCIPintervalIsEmpty(infinity, operand));
2641 
2642  intervallen = operand.sup - operand.inf;
2643  if( intervallen >= 2*M_PI )
2644  {
2645  SCIPintervalSetBounds(resultant, -1.0, 1.0);
2646  return;
2647  }
2648 
2649  modinf = fmod(operand.inf, 2*M_PI);
2650  if( modinf < 0.0 )
2651  modinf += 2*M_PI;
2652  modsup = modinf + intervallen;
2653 
2654  for( b = 0; ; ++b )
2655  {
2656  if( modinf <= sin_extremepoints[b] )
2657  {
2658  a = b;
2659  break;
2660  }
2661  }
2662  for( ; b < 4; ++b )
2663  {
2664  if( modsup <= sin_extremepoints[b] )
2665  break;
2666  }
2667 
2668  nbetween = b-a;
2669 
2670  if( nbetween > 1 )
2671  {
2672  SCIPintervalSetBounds(resultant, -1.0, 1.0);
2673  return;
2674  }
2675 
2676  finf = sin(operand.inf);
2677  fsup = sin(operand.sup);
2678 
2679  if( nbetween == 0 )
2680  {
2681  if( a & 1 ) /* next extremepoint is minimum -> decreasing -> finf < fsup */
2682  SCIPintervalSetBounds(resultant, fsup, finf);
2683  else
2684  SCIPintervalSetBounds(resultant, finf, fsup);
2685  }
2686  else /* 1 extremepoint in between */
2687  {
2688  if( a & 1 ) /* extremepoint is minimum */
2689  SCIPintervalSetBounds(resultant, -1.0, MAX(finf,fsup));
2690  else
2691  SCIPintervalSetBounds(resultant, MIN(finf,fsup), 1.0);
2692  }
2693 
2694  /* above operations did not take outward rounding into account,
2695  * so extend the computed interval slightly to increase the chance that it will contain the complete sin(operand)
2696  */
2697  if( resultant->inf > -1.0 )
2698  resultant->inf = MAX(-1.0, resultant->inf - 1e-10 * REALABS(resultant->inf)); /*lint !e666*/
2699  if( resultant->sup < 1.0 )
2700  resultant->sup = MIN( 1.0, resultant->sup + 1e-10 * REALABS(resultant->sup)); /*lint !e666*/
2701 
2702  assert(resultant->inf <= resultant->sup);
2703 }
2704 
2705 /** stores cosine value of operand in resultant
2706  * NOTE: the operations are not applied rounding-safe here
2707  */
2709  SCIP_Real infinity, /**< value for infinity */
2710  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2711  SCIP_INTERVAL operand /**< operand of operation */
2712  )
2713 {
2714  SCIP_Real intervallen;
2715  SCIP_Real modinf;
2716  SCIP_Real modsup;
2717  SCIP_Real finf;
2718  SCIP_Real fsup;
2719  int a;
2720  int b;
2721  int nbetween;
2722 
2723  assert(resultant != NULL);
2724  assert(!SCIPintervalIsEmpty(infinity, operand));
2725 
2726  intervallen = operand.sup - operand.inf;
2727  if( intervallen >= 2*M_PI )
2728  {
2729  SCIPintervalSetBounds(resultant, -1.0, 1.0);
2730  return;
2731  }
2732 
2733  modinf = fmod(operand.inf, 2*M_PI);
2734  if( modinf < 0.0 )
2735  modinf += 2*M_PI;
2736  modsup = modinf + intervallen;
2737 
2738  for( b = 0; ; ++b )
2739  {
2740  if( modinf <= cos_extremepoints[b] )
2741  {
2742  a = b;
2743  break;
2744  }
2745  }
2746  for( ; b < 3; ++b )
2747  {
2748  if( modsup <= cos_extremepoints[b] )
2749  break;
2750  }
2751 
2752  nbetween = b-a;
2753 
2754  if( nbetween > 1 )
2755  {
2756  SCIPintervalSetBounds(resultant, -1.0, 1.0);
2757  return;
2758  }
2759 
2760  finf = cos(operand.inf);
2761  fsup = cos(operand.sup);
2762 
2763  if( nbetween == 0 )
2764  {
2765  if( a & 1 ) /* next extremepoint is maximum -> increasing -> finf < fsup */
2766  SCIPintervalSetBounds(resultant, finf, fsup);
2767  else
2768  SCIPintervalSetBounds(resultant, fsup, finf);
2769  }
2770  else /* 1 extremepoint in between */
2771  {
2772  if( a & 1 ) /* extremepoint is maximum */
2773  SCIPintervalSetBounds(resultant, MIN(finf,fsup), 1.0);
2774  else
2775  SCIPintervalSetBounds(resultant, -1.0, MAX(finf,fsup));
2776  }
2777 
2778  /* above operations did not take outward rounding into account,
2779  * so extend the computed interval slightly to increase the chance that it will contain the complete cos(operand)
2780  */
2781  if( resultant->inf > -1.0 )
2782  resultant->inf = MAX(-1.0, resultant->inf - 1e-10 * REALABS(resultant->inf)); /*lint !e666*/
2783  if( resultant->sup < 1.0 )
2784  resultant->sup = MIN( 1.0, resultant->sup + 1e-10 * REALABS(resultant->sup)); /*lint !e666*/
2785 
2786  assert(resultant->inf <= resultant->sup);
2787 }
2788 
2789 /** stores sign of operand in resultant */
2791  SCIP_Real infinity, /**< value for infinity */
2792  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2793  SCIP_INTERVAL operand /**< operand of operation */
2794  )
2795 {
2796  assert(resultant != NULL);
2797  assert(!SCIPintervalIsEmpty(infinity, operand));
2798 
2799  if( operand.sup < 0.0 )
2800  {
2801  resultant->inf = -1.0;
2802  resultant->sup = -1.0;
2803  }
2804  else if( operand.inf >= 0.0 )
2805  {
2806  resultant->inf = 1.0;
2807  resultant->sup = 1.0;
2808  }
2809  else
2810  {
2811  resultant->inf = -1.0;
2812  resultant->sup = 1.0;
2813  }
2814 }
2815 
2816 /** computes exact upper bound on \f$ a x^2 + b x \f$ for x in [xlb, xub], b an interval, and a scalar
2817  *
2818  * Uses Algorithm 2.2 from Domes and Neumaier: Constraint propagation on quadratic constraints (2008) */
2820  SCIP_Real infinity, /**< value for infinity */
2821  SCIP_Real a, /**< coefficient of x^2 */
2822  SCIP_INTERVAL b_, /**< coefficient of x */
2823  SCIP_INTERVAL x /**< range of x */
2824  )
2825 {
2826  SCIP_Real b;
2827  SCIP_Real u;
2828 
2829  assert(!SCIPintervalIsEmpty(infinity, x));
2830  assert(b_.inf < infinity);
2831  assert(b_.sup > -infinity);
2832  assert( x.inf < infinity);
2833  assert( x.sup > -infinity);
2834 
2835  /* handle b*x separately */
2836  if( a == 0.0 )
2837  {
2838  if( (b_.inf <= -infinity && x.inf < 0.0 ) ||
2839  ( b_.inf < 0.0 && x.inf <= -infinity) ||
2840  ( b_.sup > 0.0 && x.sup >= infinity) ||
2841  ( b_.sup >= infinity && x.sup > 0.0 ) )
2842  {
2843  u = infinity;
2844  }
2845  else
2846  {
2847  SCIP_ROUNDMODE roundmode;
2848  SCIP_Real cand1, cand2, cand3, cand4;
2849 
2850  roundmode = intervalGetRoundingMode();
2852  cand1 = b_.inf * x.inf;
2853  cand2 = b_.inf * x.sup;
2854  cand3 = b_.sup * x.inf;
2855  cand4 = b_.sup * x.sup;
2856  u = MAX(MAX(cand1, cand2), MAX(cand3, cand4));
2857 
2858  intervalSetRoundingMode(roundmode);
2859  }
2860 
2861  return u;
2862  }
2863 
2864  if( x.sup <= 0.0 )
2865  { /* change sign of x: enclose a*x^2 + [-bub, -blb]*(-x) for (-x) in [-xub, -xlb] */
2866  u = x.sup;
2867  x.sup = -x.inf;
2868  x.inf = -u;
2869  b = -b_.inf;
2870  }
2871  else
2872  {
2873  b = b_.sup;
2874  }
2875 
2876  if( x.inf >= 0.0 )
2877  { /* upper bound for a*x^2 + b*x */
2878  SCIP_ROUNDMODE roundmode;
2879  SCIP_Real s,t;
2880 
2881  if( b >= infinity )
2882  return infinity;
2883 
2884  roundmode = intervalGetRoundingMode();
2886 
2887  u = MAX(x.inf * (a*x.inf + b), x.sup * (a*x.sup + b));
2888  s = b/2;
2889  t = s/negate(a);
2890  if( t > x.inf && negate(2*a)*x.sup > b && s*t > u )
2891  u = s*t;
2892 
2893  intervalSetRoundingMode(roundmode);
2894  return u;
2895  }
2896  else
2897  {
2898  SCIP_INTERVAL xlow = x;
2899  SCIP_Real cand1;
2900  SCIP_Real cand2;
2901  assert(x.inf < 0.0 && x.sup > 0);
2902 
2903  xlow.sup = 0; /* so xlow is lower part of interval */
2904  x.inf = 0; /* so x is upper part of interval now */
2905  cand1 = SCIPintervalQuadUpperBound(infinity, a, b_, xlow);
2906  cand2 = SCIPintervalQuadUpperBound(infinity, a, b_, x);
2907  return MAX(cand1, cand2);
2908  }
2909 }
2910 
2911 /** stores range of quadratic term in resultant
2912  *
2913  * given scalar a and intervals b and x, computes interval for \f$ a x^2 + b x \f$ */
2915  SCIP_Real infinity, /**< value for infinity */
2916  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2917  SCIP_Real sqrcoeff, /**< coefficient of x^2 */
2918  SCIP_INTERVAL lincoeff, /**< coefficient of x */
2919  SCIP_INTERVAL xrng /**< range of x */
2920  )
2921 {
2922  SCIP_Real tmp;
2923 
2924  if( SCIPintervalIsEmpty(infinity, xrng) )
2925  {
2926  SCIPintervalSetEmpty(resultant);
2927  return;
2928  }
2929  if( sqrcoeff == 0.0 )
2930  {
2931  SCIPintervalMul(infinity, resultant, lincoeff, xrng);
2932  return;
2933  }
2934 
2935  resultant->sup = SCIPintervalQuadUpperBound(infinity, sqrcoeff, lincoeff, xrng);
2936 
2937  tmp = lincoeff.inf;
2938  lincoeff.inf = -lincoeff.sup;
2939  lincoeff.sup = -tmp;
2940  resultant->inf = -SCIPintervalQuadUpperBound(infinity, -sqrcoeff, lincoeff, xrng);
2941 
2942  assert(resultant->sup >= resultant->inf);
2943 }
2944 
2945 /** computes interval with positive solutions of a quadratic equation with interval coefficients
2946  *
2947  * Given intervals a, b, and c, this function computes an interval that contains all positive solutions of \f$ a x^2 + b x \in c\f$ within xbnds.
2948  */
2950  SCIP_Real infinity, /**< value for infinity */
2951  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2952  SCIP_INTERVAL sqrcoeff, /**< coefficient of x^2 */
2953  SCIP_INTERVAL lincoeff, /**< coefficient of x */
2954  SCIP_INTERVAL rhs, /**< right hand side of equation */
2955  SCIP_INTERVAL xbnds /**< bounds on x */
2956  )
2957 {
2958  assert(resultant != NULL);
2959 
2960  /* find x>=0 s.t. a.inf x^2 + b.inf x <= c.sup -> -a.inf x^2 - b.inf x >= -c.sup */
2961  if( lincoeff.inf <= -infinity || rhs.sup >= infinity || sqrcoeff.inf <= -infinity )
2962  {
2963  resultant->inf = 0.0;
2964  resultant->sup = infinity;
2965  }
2966  else
2967  {
2968  SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(infinity, resultant, -sqrcoeff.inf, -lincoeff.inf, -rhs.sup, xbnds);
2969  SCIPdebugMessage("solve %g*x^2 + %g*x >= %g gives [%.20f, %.20f]\n", -sqrcoeff.inf, -lincoeff.inf, -rhs.sup, resultant->inf, resultant->sup);
2970  }
2971 
2972  /* find x>=0 s.t. a.sup x^2 + b.sup x >= c.inf */
2973  if( lincoeff.sup < infinity && rhs.inf > -infinity && sqrcoeff.sup < infinity ) /*lint !e644*/
2974  {
2975  SCIP_INTERVAL res2;
2976  /* coverity[uninit_use_in_call] */
2977  SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(infinity, &res2, sqrcoeff.sup, lincoeff.sup, rhs.inf, xbnds); /*lint !e644*/
2978  SCIPdebugMessage("solve %g*x^2 + %g*x >= %g gives [%.20f, %.20f]\n", sqrcoeff.sup, lincoeff.sup, rhs.inf, res2.inf, res2.sup);
2979  SCIPdebugMessage("intersection of [%.20f, %.20f] and [%.20f, %.20f]", resultant->inf, resultant->sup, res2.inf, res2.sup);
2980  /* intersect both results */
2981  SCIPintervalIntersect(resultant, *resultant, res2);
2982  SCIPdebugPrintf(" gives [%.20f, %.20f]\n", resultant->inf, resultant->sup);
2983  }
2984  /* else res2 = [0, infty] */
2985 
2986  if( resultant->inf >= infinity || resultant->sup <= -infinity )
2987  {
2988  SCIPintervalSetEmpty(resultant);
2989  }
2990 }
2991 
2992 /** computes interval with negative solutions of a quadratic equation with interval coefficients
2993  *
2994  * Given intervals a, b, and c, this function computes an interval that contains all negative solutions of \f$ a x^2 + b x \in c\f$ within xbnds.
2995  */
2997  SCIP_Real infinity, /**< value for infinity */
2998  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2999  SCIP_INTERVAL sqrcoeff, /**< coefficient of x^2 */
3000  SCIP_INTERVAL lincoeff, /**< coefficient of x */
3001  SCIP_INTERVAL rhs, /**< right hand side of equation */
3002  SCIP_INTERVAL xbnds /**< bounds on x */
3003  )
3004 {
3005  SCIP_Real tmp;
3006 
3007  /* change in variables y = -x, thus get all positive solutions of
3008  * a * y^2 + (-b) * y in c with -xbnds as bounds on y
3009  */
3010 
3011  tmp = lincoeff.inf;
3012  lincoeff.inf = -lincoeff.sup;
3013  lincoeff.sup = -tmp;
3014 
3015  tmp = xbnds.inf;
3016  xbnds.inf = -xbnds.sup;
3017  xbnds.sup = -tmp;
3018 
3019  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, resultant, sqrcoeff, lincoeff, rhs, xbnds);
3020 
3021  tmp = resultant->inf;
3022  resultant->inf = -resultant->sup;
3023  resultant->sup = -tmp;
3024 }
3025 
3026 
3027 /** computes positive solutions of a quadratic equation with scalar coefficients
3028  *
3029  * Given scalar a, b, and c, this function computes an interval that contains all positive solutions of \f$ a x^2 + b x \geq c\f$ within xbnds.
3030  * Implements Algorithm 3.2 from Domes and Neumaier: Constraint propagation on quadratic constraints (2008).
3031  */
3033  SCIP_Real infinity, /**< value for infinity */
3034  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
3035  SCIP_Real sqrcoeff, /**< coefficient of x^2 */
3036  SCIP_Real lincoeff, /**< coefficient of x */
3037  SCIP_Real rhs, /**< right hand side of equation */
3038  SCIP_INTERVAL xbnds /**< bounds on x */
3039  )
3040 {
3041  SCIP_ROUNDMODE roundmode;
3042  SCIP_Real b;
3043  SCIP_Real delta;
3044  SCIP_Real z;
3045 
3046  assert(resultant != NULL);
3047  assert(sqrcoeff < infinity);
3048  assert(sqrcoeff > -infinity);
3049 
3050  if( sqrcoeff == 0.0 )
3051  {
3052  /* special handling for linear b * x >= c
3053  *
3054  * The non-negative solutions here are:
3055  * b < 0, c <= 0 : [0, c/b]
3056  * b <= 0, c > 0 : empty
3057  * b > 0, c > 0 : [c/b, infty]
3058  * b >= 0, c <= 0 : [0, infty]
3059  *
3060  * The same should have been computed below, but without the sqrcoeff, terms simplify (thus, also less rounding).
3061  */
3062 
3063  if( lincoeff <= 0.0 && rhs > 0.0 )
3064  {
3065  SCIPintervalSetEmpty(resultant);
3066  return;
3067  }
3068 
3069  if( lincoeff >= 0.0 && rhs <= 0.0 )
3070  {
3071  /* [0,infty] cap xbnds */
3072  resultant->inf = MAX(0.0, xbnds.inf);
3073  resultant->sup = xbnds.sup;
3074  return;
3075  }
3076 
3077  roundmode = intervalGetRoundingMode();
3078 
3079  if( lincoeff < 0.0 && rhs <= 0.0 )
3080  {
3081  /* [0,c/b] cap xbnds */
3082  resultant->inf = MAX(0.0, xbnds.inf);
3083 
3085  resultant->sup = rhs / lincoeff;
3086  if( xbnds.sup < resultant->sup )
3087  resultant->sup = xbnds.sup;
3088  }
3089  else
3090  {
3091  assert(lincoeff > 0.0);
3092  assert(rhs > 0.0);
3093 
3094  /* [c/b, infty] cap xbnds */
3095 
3097  resultant->inf = rhs / lincoeff;
3098  if( resultant->inf < xbnds.inf )
3099  resultant->inf = xbnds.inf;
3100 
3101  resultant->sup = xbnds.sup;
3102  }
3103 
3104  intervalSetRoundingMode(roundmode);
3105 
3106  return;
3107  }
3108 
3109  resultant->inf = 0.0;
3110  resultant->sup = infinity;
3111 
3112  roundmode = intervalGetRoundingMode();
3113 
3114  /* this should actually be round_upwards, but unless lincoeff is min_double,
3115  * there shouldn't be any rounding happening when dividing by 2, i.e., shifting exponent,
3116  * so it is ok to not change the rounding mode here
3117  */
3118  b = lincoeff / 2.0;
3119 
3120  if( lincoeff >= 0.0 )
3121  { /* b >= 0.0 */
3122  if( rhs > 0.0 )
3123  { /* b >= 0.0 and c > 0.0 */
3125  delta = b*b + sqrcoeff*rhs;
3126  if( delta < 0.0 )
3127  {
3128  SCIPintervalSetEmpty(resultant);
3129  }
3130  else
3131  {
3133  z = SCIPnextafter(sqrt(delta), SCIP_REAL_MAX);
3135  z += b;
3136  resultant->inf = negate(negate(rhs)/z);
3137  if( sqrcoeff < 0.0 )
3138  resultant->sup = z / negate(sqrcoeff);
3139  }
3140  }
3141  else
3142  { /* b >= 0.0 and c <= 0.0 */
3143  if( sqrcoeff < 0.0 )
3144  {
3146  delta = b*b + sqrcoeff*rhs;
3148  z = SCIPnextafter(sqrt(delta), SCIP_REAL_MAX);
3150  z += b;
3151  resultant->sup = z / negate(sqrcoeff);
3152  }
3153  }
3154  }
3155  else
3156  { /* b < 0.0 */
3157  if( rhs > 0.0 )
3158  { /* b < 0.0 and c > 0.0 */
3159  if( sqrcoeff > 0.0 )
3160  {
3162  delta = b*b + sqrcoeff*rhs;
3164  z = SCIPnextafter(sqrt(delta), SCIP_REAL_MIN);
3166  z += negate(b);
3167  resultant->inf = z / sqrcoeff;
3168  }
3169  else
3170  {
3171  SCIPintervalSetEmpty(resultant);
3172  }
3173  }
3174  else
3175  { /* b < 0.0 and c <= 0.0 */
3177  delta = b*b + sqrcoeff * rhs;
3178  if( delta >= 0.0 )
3179  {
3180  /* let resultant = [0,-c/z] for now */
3182  z = SCIPnextafter(sqrt(delta), SCIP_REAL_MIN);
3183  /* continue with downward rounding, because we want z (>= 0) to be small,
3184  * because -rhs/z needs to be large (-rhs >= 0)
3185  */
3187  z += negate(b);
3188  /* also now compute rhs/z with downward rounding, so that -(rhs/z) becomes large */
3189  resultant->sup = negate(rhs/z);
3190 
3191  if( sqrcoeff > 0.0 )
3192  {
3193  /* for a > 0, the result is [0,-c/z] \vee [z/a,infinity]
3194  * currently, resultant = [0,-c/z]
3195  */
3196  SCIP_Real zdiva;
3197 
3198  zdiva = z/sqrcoeff;
3199 
3200  if( xbnds.sup < zdiva )
3201  {
3202  /* after intersecting with xbnds, result is [0,-c/z], so we are done */
3203  }
3204  else if( xbnds.inf > resultant->sup )
3205  {
3206  /* after intersecting with xbnds, result is [z/a,infinity] */
3207  resultant->inf = zdiva;
3208  resultant->sup = infinity;
3209  }
3210  else
3211  {
3212  /* after intersecting with xbnds we can neither exclude [0,-c/z] nor [z/a,infinity],
3213  * so put resultant = [0,infinity] (intersection with xbnds happens below)
3214  * @todo we could create a hole here
3215  */
3216  resultant->sup = infinity;
3217  }
3218  }
3219  else
3220  {
3221  /* for a < 0, the result is [0,-c/z], so we are done */
3222  }
3223  }
3224  }
3225  }
3226 
3227  SCIPintervalIntersect(resultant, *resultant, xbnds);
3228 
3229  intervalSetRoundingMode(roundmode);
3230 }
3231 
3232 /** solves a quadratic equation with interval coefficients
3233  *
3234  * Given intervals a, b and c, this function computes an interval that contains all solutions of \f$ a x^2 + b x \in c\f$ within xbnds
3235  */
3237  SCIP_Real infinity, /**< value for infinity */
3238  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
3239  SCIP_INTERVAL sqrcoeff, /**< coefficient of x^2 */
3240  SCIP_INTERVAL lincoeff, /**< coefficient of x */
3241  SCIP_INTERVAL rhs, /**< right hand side of equation */
3242  SCIP_INTERVAL xbnds /**< bounds on x */
3243  )
3244 {
3245  SCIP_INTERVAL xpos;
3246  SCIP_INTERVAL xneg;
3247 
3248  assert(resultant != NULL);
3249  assert(!SCIPintervalIsEmpty(infinity, sqrcoeff));
3250  assert(!SCIPintervalIsEmpty(infinity, lincoeff));
3251  assert(!SCIPintervalIsEmpty(infinity, rhs));
3252 
3253  /* special handling for lincoeff * x = rhs without 0 in lincoeff
3254  * rhs/lincoeff gives a good interval that we just have to intersect with xbnds
3255  * the code below would also work, but uses many more case distinctions to get to a result that should be the same (though epsilon differences can sometimes be observed)
3256  */
3257  if( sqrcoeff.inf == 0.0 && sqrcoeff.sup == 0.0 && (lincoeff.inf > 0.0 || lincoeff.sup < 0.0) )
3258  {
3259  SCIPintervalDiv(infinity, resultant, rhs, lincoeff);
3260  SCIPintervalIntersect(resultant, *resultant, xbnds);
3261  SCIPdebugMessage("solving [%g,%g]*x = [%g,%g] for x in [%g,%g] gives [%g,%g]\n", lincoeff.inf, lincoeff.sup, rhs.inf, rhs.sup, xbnds.inf, xbnds.sup, resultant->inf, resultant->sup);
3262  return;
3263  }
3264 
3265  SCIPdebugMessage("solving [%g,%g]*x^2 + [%g,%g]*x = [%g,%g] for x in [%g,%g]\n", sqrcoeff.inf, sqrcoeff.sup, lincoeff.inf, lincoeff.sup, rhs.inf, rhs.sup, xbnds.inf, xbnds.sup);
3266 
3267  /* find all x>=0 such that a*x^2+b*x = c */
3268  if( xbnds.sup >= 0 )
3269  {
3270  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, &xpos, sqrcoeff, lincoeff, rhs, xbnds);
3271  SCIPdebugMessage(" solutions of [%g,%g]*x^2 + [%g,%g]*x in [%g,%g] for x in [%g,%g] are [%.15g,%.15g]\n",
3272  sqrcoeff.inf, sqrcoeff.sup, lincoeff.inf, lincoeff.sup, rhs.inf, rhs.sup, MAX(xbnds.inf, 0.0), xbnds.sup, xpos.inf, xpos.sup);
3273  }
3274  else
3275  {
3276  SCIPintervalSetEmpty(&xpos);
3277  }
3278 
3279  /* find all x<=0 such that a*x^2-b*x = c */
3280  if( xbnds.inf <= 0.0 )
3281  {
3282  SCIPintervalSolveUnivariateQuadExpressionNegative(infinity, &xneg, sqrcoeff, lincoeff, rhs, xbnds);
3283  SCIPdebugMessage(" solutions of [%g,%g]*x^2 + [%g,%g]*x in [%g,%g] for x in [%g,%g] are [%g,%g]\n",
3284  sqrcoeff.inf, sqrcoeff.sup, lincoeff.inf, lincoeff.sup, rhs.inf, rhs.sup, xbnds.inf, MIN(xbnds.sup, 0.0), xneg.inf, xneg.sup);
3285  }
3286  else
3287  {
3288  SCIPintervalSetEmpty(&xneg);
3289  }
3290 
3291  SCIPintervalUnify(resultant, xpos, xneg);
3292  SCIPdebugMessage(" unify gives [%g,%g]\n", SCIPintervalGetInf(*resultant), SCIPintervalGetSup(*resultant));
3293 }
3294 
3295 /** stores range of bivariate quadratic term in resultant
3296  * given scalars ax, ay, axy, bx, and by and intervals for x and y, computes interval for \f$ ax x^2 + ay y^2 + axy x y + bx x + by y \f$
3297  * NOTE: the operations are not applied rounding-safe here
3298  */
3300  SCIP_Real infinity, /**< value for infinity in interval arithmetics */
3301  SCIP_INTERVAL* resultant, /**< buffer where to store result of operation */
3302  SCIP_Real ax, /**< square coefficient of x */
3303  SCIP_Real ay, /**< square coefficient of y */
3304  SCIP_Real axy, /**< bilinear coefficients */
3305  SCIP_Real bx, /**< linear coefficient of x */
3306  SCIP_Real by, /**< linear coefficient of y */
3307  SCIP_INTERVAL xbnds, /**< bounds on x */
3308  SCIP_INTERVAL ybnds /**< bounds on y */
3309  )
3310 {
3311  /* we use double double precision and finally widen the computed range by 1e-8% to compensate for not computing rounding-safe here */
3312  SCIP_Real minval;
3313  SCIP_Real maxval;
3314  SCIP_Real val;
3315  SCIP_Real x;
3316  SCIP_Real y;
3317  SCIP_Real denom;
3318 
3319  assert(resultant != NULL);
3320  assert(xbnds.inf <= xbnds.sup);
3321  assert(ybnds.inf <= ybnds.sup);
3322 
3323  /* if we are separable, then fall back to use SCIPintervalQuad two times and add */
3324  if( axy == 0.0 )
3325  {
3326  SCIP_INTERVAL tmp;
3327 
3328  SCIPintervalSet(&tmp, bx);
3329  SCIPintervalQuad(infinity, resultant, ax, tmp, xbnds);
3330 
3331  SCIPintervalSet(&tmp, by);
3332  SCIPintervalQuad(infinity, &tmp, ay, tmp, ybnds);
3333 
3334  SCIPintervalAdd(infinity, resultant, *resultant, tmp);
3335 
3336  return;
3337  }
3338 
3339  SCIPintervalSet(resultant, 0.0);
3340 
3341  minval = infinity;
3342  maxval = -infinity;
3343 
3344  /* check minima/maxima of expression */
3345  denom = 4.0 * ax * ay - axy * axy;
3346  if( REALABS(denom) > 1e-9 )
3347  {
3348  x = (axy * by - 2.0 * ay * bx) / denom;
3349  y = (axy * bx - 2.0 * ax * by) / denom;
3350  if( xbnds.inf <= x && x <= xbnds.sup && ybnds.inf <= y && y <= ybnds.sup )
3351  {
3352  val = (axy * bx * by - ay * bx * bx - ax * by * by) / denom;
3353  minval = MIN(val, minval);
3354  maxval = MAX(val, maxval);
3355  }
3356  }
3357  else if( REALABS(2.0 * ay * bx - axy * by) <= 1e-9 )
3358  {
3359  /* The whole line (x, -bx/axy - (axy/2ay) x) defines an extreme point with value -ay bx^2 / axy^2
3360  * If x is unbounded, then there is an (x,y) with y in ybnds where the extreme value is assumed.
3361  * If x is bounded on at least one side, then we can rely that the checks below for x at one of its bounds will check this extreme point.
3362  */
3363  if( xbnds.inf <= -infinity && xbnds.sup >= infinity )
3364  {
3365  val = -ay * bx * bx / (axy * axy);
3366  minval = MIN(val, minval);
3367  maxval = MAX(val, maxval);
3368  }
3369  }
3370 
3371  /* check boundary of box xbnds x ybnds */
3372 
3373  if( xbnds.inf <= -infinity )
3374  {
3375  /* check value for x -> -infinity */
3376  if( ax > 0.0 )
3377  maxval = infinity;
3378  else if( ax < 0.0 )
3379  minval = -infinity;
3380  else if( ax == 0.0 )
3381  {
3382  /* bivar(x,y) tends to -(bx+axy y) * infinity */
3383 
3384  if( ybnds.inf <= -infinity )
3385  val = (axy < 0.0 ? -infinity : infinity);
3386  else if( bx + axy * ybnds.inf < 0.0 )
3387  val = infinity;
3388  else
3389  val = -infinity;
3390  minval = MIN(val, minval);
3391  maxval = MAX(val, maxval);
3392 
3393  if( ybnds.sup >= infinity )
3394  val = (axy < 0.0 ? infinity : -infinity);
3395  else if( bx + axy * ybnds.sup < 0.0 )
3396  val = infinity;
3397  else
3398  val = -infinity;
3399  minval = MIN(val, minval);
3400  maxval = MAX(val, maxval);
3401  }
3402  }
3403  else
3404  {
3405  /* get range of bivar(xbnds.inf, y) for y in ybnds */
3406  SCIP_INTERVAL tmp;
3407  SCIP_INTERVAL ycoef;
3408 
3409  SCIPintervalSet(&ycoef, axy * xbnds.inf + by);
3410  SCIPintervalQuad(infinity, &tmp, ay, ycoef, ybnds);
3411  SCIPintervalAddScalar(infinity, &tmp, tmp, (SCIP_Real)(ax * xbnds.inf * xbnds.inf + bx * xbnds.inf));
3412  minval = MIN(tmp.inf, minval);
3413  maxval = MAX(tmp.sup, maxval);
3414  }
3415 
3416  if( xbnds.sup >= infinity )
3417  {
3418  /* check value for x -> infinity */
3419  if( ax > 0.0 )
3420  maxval = infinity;
3421  else if( ax < 0.0 )
3422  minval = -infinity;
3423  else if( ax == 0.0 )
3424  {
3425  /* bivar(x,y) tends to (bx+axy y) * infinity */
3426 
3427  if( ybnds.inf <= -infinity )
3428  val = (axy > 0.0 ? -infinity : infinity);
3429  else if( bx + axy * ybnds.inf > 0.0 )
3430  val = infinity;
3431  else
3432  val = -infinity;
3433  minval = MIN(val, minval);
3434  maxval = MAX(val, maxval);
3435 
3436  if( ybnds.sup >= infinity )
3437  val = (axy > 0.0 ? infinity : -infinity);
3438  else if( bx + axy * ybnds.sup > 0.0 )
3439  val = infinity;
3440  else
3441  val = -infinity;
3442  minval = MIN(val, minval);
3443  maxval = MAX(val, maxval);
3444  }
3445  }
3446  else
3447  {
3448  /* get range of bivar(xbnds.sup, y) for y in ybnds */
3449  SCIP_INTERVAL tmp;
3450  SCIP_INTERVAL ycoef;
3451 
3452  SCIPintervalSet(&ycoef, axy * xbnds.sup + by);
3453  SCIPintervalQuad(infinity, &tmp, ay, ycoef, ybnds);
3454  SCIPintervalAddScalar(infinity, &tmp, tmp, (SCIP_Real)(ax * xbnds.sup * xbnds.sup + bx * xbnds.sup));
3455  minval = MIN(tmp.inf, minval);
3456  maxval = MAX(tmp.sup, maxval);
3457  }
3458 
3459  if( ybnds.inf <= -infinity )
3460  {
3461  /* check value for y -> -infinity */
3462  if( ay > 0.0 )
3463  maxval = infinity;
3464  else if( ay < 0.0 )
3465  minval = -infinity;
3466  else if( ay == 0.0 )
3467  {
3468  /* bivar(x,y) tends to -(by+axy x) * infinity */
3469 
3470  if( xbnds.inf <= -infinity )
3471  val = (axy < 0.0 ? -infinity : infinity);
3472  else if( by + axy * xbnds.inf < 0.0 )
3473  val = infinity;
3474  else
3475  val = -infinity;
3476  minval = MIN(val, minval);
3477  maxval = MAX(val, maxval);
3478 
3479  if( xbnds.sup >= infinity )
3480  val = (axy < 0.0 ? infinity : -infinity);
3481  else if( by + axy * xbnds.sup < 0.0 )
3482  val = infinity;
3483  else
3484  val = -infinity;
3485  minval = MIN(val, minval);
3486  maxval = MAX(val, maxval);
3487  }
3488  }
3489  else
3490  {
3491  /* get range of bivar(x, ybnds.inf) for x in xbnds */
3492  SCIP_INTERVAL tmp;
3493  SCIP_INTERVAL xcoef;
3494 
3495  SCIPintervalSet(&xcoef, axy * ybnds.inf + bx);
3496  SCIPintervalQuad(infinity, &tmp, ax, xcoef, xbnds);
3497  SCIPintervalAddScalar(infinity, &tmp, tmp, ay * ybnds.inf * ybnds.inf + by * ybnds.inf);
3498  minval = MIN(tmp.inf, minval);
3499  maxval = MAX(tmp.sup, maxval);
3500  }
3501 
3502  if( ybnds.sup >= infinity )
3503  {
3504  /* check value for y -> infinity */
3505  if( ay > 0.0 )
3506  maxval = infinity;
3507  else if( ay < 0.0 )
3508  minval = -infinity;
3509  else if( ay == 0.0 )
3510  {
3511  /* bivar(x,y) tends to (by+axy x) * infinity */
3512 
3513  if( xbnds.inf <= -infinity )
3514  val = (axy > 0.0 ? -infinity : infinity);
3515  else if( by + axy * xbnds.inf > 0.0 )
3516  val = infinity;
3517  else
3518  val = -infinity;
3519  minval = MIN(val, minval);
3520  maxval = MAX(val, maxval);
3521 
3522  if( xbnds.sup >= infinity )
3523  val = (axy > 0.0 ? infinity : -infinity);
3524  else if( by + axy * xbnds.sup > 0.0 )
3525  val = infinity;
3526  else
3527  val = -infinity;
3528  minval = MIN(val, minval);
3529  maxval = MAX(val, maxval);
3530  }
3531  }
3532  else
3533  {
3534  /* get range of bivar(x, ybnds.sup) for x in xbnds */
3535  SCIP_INTERVAL tmp;
3536  SCIP_INTERVAL xcoef;
3537 
3538  SCIPintervalSet(&xcoef, axy * ybnds.sup + bx);
3539  SCIPintervalQuad(infinity, &tmp, ax, xcoef, xbnds);
3540  SCIPintervalAddScalar(infinity, &tmp, tmp, (SCIP_Real)(ay * ybnds.sup * ybnds.sup + by * ybnds.sup));
3541  minval = MIN(tmp.inf, minval);
3542  maxval = MAX(tmp.sup, maxval);
3543  }
3544 
3545  minval -= 1e-10 * REALABS(minval);
3546  maxval += 1e-10 * REALABS(maxval);
3547  SCIPintervalSetBounds(resultant, (SCIP_Real)minval, (SCIP_Real)maxval);
3548 
3549  SCIPdebugMessage("range for %gx^2 + %gy^2 + %gxy + %gx + %gy = [%g, %g] for x = [%g, %g], y=[%g, %g]\n",
3550  ax, ay, axy, bx, by, minval, maxval, xbnds.inf, xbnds.sup, ybnds.inf, ybnds.sup);
3551 }
3552 
3553 /** solves a bivariate quadratic equation for the first variable
3554  * given scalars ax, ay, axy, bx and by, and intervals for x, y, and rhs,
3555  * computes \f$ \{ x \in \mathbf{x} : \exists y \in \mathbf{y} : a_x x^2 + a_y y^2 + a_{xy} x y + b_x x + b_y y \in \mathbf{\mbox{rhs}} \} \f$
3556  * NOTE: the operations are not applied rounding-safe here
3557  */
3559  SCIP_Real infinity, /**< value for infinity in interval arithmetics */
3560  SCIP_INTERVAL* resultant, /**< buffer where to store result of operation */
3561  SCIP_Real ax, /**< square coefficient of x */
3562  SCIP_Real ay, /**< square coefficient of y */
3563  SCIP_Real axy, /**< bilinear coefficients */
3564  SCIP_Real bx, /**< linear coefficient of x */
3565  SCIP_Real by, /**< linear coefficient of y */
3566  SCIP_INTERVAL rhs, /**< right-hand-side of equation */
3567  SCIP_INTERVAL xbnds, /**< bounds on x */
3568  SCIP_INTERVAL ybnds /**< bounds on y */
3569  )
3570 {
3571  /* we use double double precision and finally widen the computed range by 1e-8% to compensate for not computing rounding-safe here */
3572  SCIP_Real val;
3573 
3574  assert(resultant != NULL);
3575 
3576  if( axy == 0.0 )
3577  {
3578  /* if axy == 0, fall back to SCIPintervalSolveUnivariateQuadExpression */
3579  SCIP_INTERVAL ytermrng;
3580  SCIP_INTERVAL sqrcoef;
3581  SCIP_INTERVAL lincoef;
3582  SCIP_INTERVAL pos;
3583  SCIP_INTERVAL neg;
3584 
3585  SCIPintervalSet(&lincoef, by);
3586  SCIPintervalQuad(infinity, &ytermrng, ay, lincoef, ybnds);
3587  SCIPintervalSub(infinity, &rhs, rhs, ytermrng);
3588 
3589  SCIPintervalSet(&sqrcoef, ax);
3590 
3591  /* get positive solutions, if of interest */
3592  if( xbnds.sup >= 0.0 )
3593  {
3594  SCIPintervalSet(&lincoef, bx);
3595  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, &pos, sqrcoef, lincoef, rhs, xbnds);
3596  }
3597  else
3598  SCIPintervalSetEmpty(&pos);
3599 
3600  /* get negative solutions, if of interest */
3601  if( xbnds.inf < 0.0 )
3602  {
3603  SCIP_INTERVAL xbndsneg;
3604  SCIPintervalSet(&lincoef, -bx);
3605  SCIPintervalSetBounds(&xbndsneg, -xbnds.sup, -xbnds.inf);
3606  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, &neg, sqrcoef, lincoef, rhs, xbndsneg);
3607  if( !SCIPintervalIsEmpty(infinity, neg) )
3608  SCIPintervalSetBounds(&neg, -neg.sup, -neg.inf);
3609  }
3610  else
3611  SCIPintervalSetEmpty(&neg);
3612 
3613  SCIPintervalUnify(resultant, pos, neg);
3614 
3615  return;
3616  }
3617 
3618  if( ybnds.inf <= -infinity || ybnds.sup >= infinity )
3619  {
3620  /* the code below is buggy if y is unbounded, see #2250
3621  * fall back to univariate case by solving a_x x^2 + b_x x + a_y y^2 + (a_xy xbnds + b_y) y in rhs
3622  */
3623  SCIP_INTERVAL ax_;
3624  SCIP_INTERVAL bx_;
3625  SCIP_INTERVAL ycoef;
3626  SCIP_INTERVAL ytermbounds;
3627 
3628  *resultant = xbnds;
3629 
3630  /* nothing we can do here if x is unbounded (we have a_xy != 0 here) */
3631  if( xbnds.inf <= -infinity && xbnds.sup >= infinity )
3632  return;
3633 
3634  /* ycoef = axy xbnds + by */
3635  SCIPintervalMulScalar(infinity, &ycoef, xbnds, axy);
3636  SCIPintervalAddScalar(infinity, &ycoef, ycoef, by);
3637 
3638  /* get bounds on ay y^2 + (axy xbnds + by) y */
3639  SCIPintervalQuad(infinity, &ytermbounds, ay, ycoef, ybnds);
3640 
3641  /* now solve ax x^2 + bx x in rhs - ytermbounds */
3642  SCIPintervalSet(&ax_, ax);
3643  SCIPintervalSet(&bx_, bx);
3644  SCIPintervalSub(infinity, &rhs, rhs, ytermbounds);
3645  SCIPintervalSolveUnivariateQuadExpression(infinity, resultant, ax_, bx_, rhs, xbnds);
3646 
3647  return;
3648  }
3649 
3650  if( ax < 0.0 )
3651  {
3652  SCIP_Real tmp;
3653  tmp = rhs.inf;
3654  rhs.inf = -rhs.sup;
3655  rhs.sup = -tmp;
3656 
3657  SCIPintervalSolveBivariateQuadExpressionAllScalar(infinity, resultant, -ax, -ay, -axy, -bx, -by, rhs, xbnds, ybnds);
3658  return;
3659  }
3660  assert(ax >= 0.0);
3661 
3662  *resultant = xbnds;
3663 
3664  if( ax > 0.0 )
3665  {
3666  SCIP_Real sqrtax;
3667  SCIP_Real minvalleft;
3668  SCIP_Real maxvalleft;
3669  SCIP_Real minvalright;
3670  SCIP_Real maxvalright;
3671  SCIP_Real ymin;
3672  SCIP_Real rcoef_y;
3673  SCIP_Real rcoef_yy;
3674  SCIP_Real rcoef_const;
3675 
3676  sqrtax = sqrt(ax);
3677 
3678  /* rewrite equation as (sqrt(ax)x + b(y))^2 \in r(rhs,y), where
3679  * b(y) = (bx + axy y)/(2sqrt(ax)), r(rhs,y) = rhs - ay y^2 - by y + b(y)^2
3680  *
3681  * -> r(rhs,y) = bx^2/(4ax) + rhs + (axy bx/(2ax) - by)*y + (axy^2/(4ax) - ay)*y^2
3682  */
3683  rcoef_y = axy * bx / (2.0*ax) - by;
3684  rcoef_yy = axy * axy / (4.0*ax) - ay;
3685  rcoef_const = bx * bx / (4.0*ax);
3686 
3687 #define CALCB(y) ((bx + axy * (y)) / (2.0 * sqrtax))
3688 #define CALCR(c,y) (rcoef_const + (c) + (rcoef_y + rcoef_yy * (y)) * (y))
3689 
3690  /* check whether r(rhs,y) is always negative */
3691  if( rhs.sup < infinity )
3692  {
3693  SCIP_INTERVAL ycoef;
3694  SCIP_Real ub;
3695 
3696  SCIPintervalSet(&ycoef, (SCIP_Real)rcoef_y);
3697  ub = (SCIP_Real)(SCIPintervalQuadUpperBound(infinity, (SCIP_Real)rcoef_yy, ycoef, ybnds) + rhs.sup + rcoef_const);
3698 
3699  if( EPSN(ub, 1e-9) )
3700  {
3701  SCIPintervalSetEmpty(resultant);
3702  return;
3703  }
3704  else if( ub < 0.0 )
3705  {
3706  /* it looks like there will be no solution (rhs < 0), but we are very close and above operations did not take care of careful rounding
3707  * thus, we relax rhs a be feasible a bit (-ub would be sufficient, but that would put us exactly onto the boundary)
3708  * see also #1861
3709  */
3710  rhs.sup += -2.0*ub;
3711  }
3712  }
3713 
3714  /* we have sqrt(ax)x \in (-sqrt(r(rhs,y))-b(y)) \cup (sqrt(r(rhs,y))-b(y))
3715  * compute minima and maxima of both functions such that
3716  *
3717  * [minvalleft, maxvalleft ] = -sqrt(r(rhs,y))-b(y)
3718  * [minvalright, maxvalright] = sqrt(r(rhs,y))-b(y)
3719  */
3720 
3721  minvalleft = infinity;
3722  maxvalleft = -infinity;
3723  minvalright = infinity;
3724  maxvalright = -infinity;
3725 
3726  if( rhs.sup >= infinity )
3727  {
3728  /* we can't do much if rhs.sup is infinite
3729  * but we may do a bit of xbnds isn't too huge and rhs.inf > -infinity
3730  */
3731  minvalleft = -infinity;
3732  maxvalright = infinity;
3733  }
3734 
3735  /* evaluate at lower bound of y, as long as r(rhs,ylb) > 0 */
3736  if( ybnds.inf <= -infinity )
3737  {
3738  /* check limit of +/-sqrt(r(rhs,y))-b(y) for y -> -infty */
3739  if( !EPSZ(ay, 1e-9) && axy * axy >= 4.0 * ax * ay )
3740  {
3741  if( axy < 0.0 )
3742  {
3743  minvalleft = -infinity;
3744 
3745  if( ay > 0.0 )
3746  minvalright = -infinity;
3747  else
3748  maxvalright = infinity;
3749  }
3750  else
3751  {
3752  maxvalright = infinity;
3753 
3754  if( ay > 0.0 )
3755  maxvalleft = infinity;
3756  else
3757  minvalleft = -infinity;
3758  }
3759  }
3760  else if( !EPSZ(ay, 1e-9) )
3761  {
3762  /* here axy * axy < 4 * ax * ay, so need to check for zeros of r(rhs,y), which is done below */
3763  }
3764  else
3765  {
3766  /* here ay = 0.0, which gives a limit of -by/2 for -sqrt(r(rhs,y))-b(y) */
3767  minvalleft = -by / 2.0;
3768  maxvalleft = -by / 2.0;
3769  /* here ay = 0.0, which gives a limit of +infinity for sqrt(r(rhs,y))-b(y) */
3770  maxvalright = infinity;
3771  }
3772  }
3773  else
3774  {
3775  SCIP_Real b;
3776  SCIP_Real c;
3777 
3778  b = CALCB(ybnds.inf);
3779 
3780  if( rhs.sup < infinity )
3781  {
3782  c = CALCR(rhs.sup, ybnds.inf);
3783 
3784  if( c > 0.0 )
3785  {
3786  SCIP_Real sqrtc;
3787 
3788  sqrtc = sqrt(c);
3789  minvalleft = MIN(-sqrtc - b, minvalleft);
3790  maxvalright = MAX( sqrtc - b, maxvalright);
3791  }
3792  }
3793 
3794  if( rhs.inf > -infinity )
3795  {
3796  c = CALCR(rhs.inf, ybnds.inf);
3797 
3798  if( c > 0.0 )
3799  {
3800  SCIP_Real sqrtc;
3801 
3802  sqrtc = sqrt(c);
3803  maxvalleft = MAX(-sqrtc - b, maxvalleft);
3804  minvalright = MIN( sqrtc - b, minvalright);
3805  }
3806  }
3807  }
3808 
3809  /* evaluate at upper bound of y, as long as r(rhs, yub) > 0 */
3810  if( ybnds.sup >= infinity )
3811  {
3812  /* check limit of +/-sqrt(r(rhs,y))-b(y) for y -> +infty */
3813  if( !EPSZ(ay, 1e-9) && axy * axy >= 4.0 * ax * ay )
3814  {
3815  if( axy > 0.0 )
3816  {
3817  minvalleft = -infinity;
3818 
3819  if( ay > 0.0 )
3820  minvalright = -infinity;
3821  else
3822  maxvalright = infinity;
3823  }
3824  else
3825  {
3826  maxvalright = infinity;
3827 
3828  if( ay > 0.0 )
3829  maxvalleft = infinity;
3830  else
3831  minvalleft = -infinity;
3832  }
3833  }
3834  else if( !EPSZ(ay, 1e-9) )
3835  {
3836  /* here axy * axy < 4 * ax * ay, so need to check for zeros of r(rhs,y), which will happen below */
3837  }
3838  else
3839  {
3840  /* here ay = 0.0, which gives a limit of -infinity for -sqrt(r(rhs,y))-b(y) */
3841  minvalleft = -infinity;
3842  /* here ay = 0.0, which gives a limit of -by/2 for sqrt(r(rhs,y))-b(y) */
3843  minvalright = MIN(minvalright, -by / 2.0);
3844  maxvalright = MAX(maxvalright, -by / 2.0);
3845  }
3846  }
3847  else
3848  {
3849  SCIP_Real b;
3850  SCIP_Real c;
3851 
3852  b = CALCB(ybnds.sup);
3853 
3854  if( rhs.sup < infinity )
3855  {
3856  c = CALCR(rhs.sup, ybnds.sup);
3857 
3858  if( c > 0.0 )
3859  {
3860  SCIP_Real sqrtc;
3861 
3862  sqrtc = sqrt(c);
3863  minvalleft = MIN(-sqrtc - b, minvalleft);
3864  maxvalright = MAX( sqrtc - b, maxvalright);
3865  }
3866  }
3867 
3868  if( rhs.inf > -infinity )
3869  {
3870  c = CALCR(rhs.inf, ybnds.sup);
3871 
3872  if( c > 0.0 )
3873  {
3874  SCIP_Real sqrtc;
3875 
3876  sqrtc = sqrt(c);
3877  maxvalleft = MAX(-sqrtc - b, maxvalleft);
3878  minvalright = MIN( sqrtc - b, minvalright);
3879  }
3880  }
3881  }
3882 
3883  /* evaluate at ymin = y_{_,+}, if inside ybnds
3884  * if ay = 0 or 2ay*bx == axy*by, then there is no ymin */
3885  if( !EPSZ(ay, 1e-9) )
3886  {
3887  if( REALABS(axy*axy - 4.0*ax*ay) > 1e-9 )
3888  {
3889  SCIP_Real sqrtterm;
3890 
3891  if( rhs.sup < infinity )
3892  {
3893  sqrtterm = axy * axy * ay * (ay * bx * bx - axy * bx * by + ax * by * by - axy * axy * rhs.sup + 4.0 * ax * ay * rhs.sup);
3894  if( !EPSN(sqrtterm, 1e-9) )
3895  {
3896  sqrtterm = sqrt(MAX(sqrtterm, 0.0));
3897  /* check first candidate for extreme points of +/-sqrt(rhs(r,y))-b(y) */
3898  ymin = axy * ay * bx - 2.0 * ax * ay * by - sqrtterm;
3899  ymin /= ay;
3900  ymin /= 4.0 * ax * ay - axy * axy;
3901 
3902  if( ymin > ybnds.inf && ymin < ybnds.sup )
3903  {
3904  SCIP_Real b;
3905  SCIP_Real c;
3906 
3907  b = CALCB(ymin);
3908  c = CALCR(rhs.sup, ymin);
3909 
3910  if( c > 0.0 )
3911  {
3912  SCIP_Real sqrtc;
3913 
3914  sqrtc = sqrt(c);
3915  minvalleft = MIN(-sqrtc - b, minvalleft);
3916  maxvalright = MAX( sqrtc - b, maxvalright);
3917  }
3918  }
3919 
3920  /* check second candidate for extreme points of +/-sqrt(rhs(r,y))-b(y) */
3921  ymin = axy * ay * bx - 2.0 * ax * ay * by + sqrtterm;
3922  ymin /= ay;
3923  ymin /= 4.0 * ax * ay - axy * axy;
3924 
3925  if( ymin > ybnds.inf && ymin < ybnds.sup )
3926  {
3927  SCIP_Real b;
3928  SCIP_Real c;
3929 
3930  b = CALCB(ymin);
3931  c = CALCR(rhs.sup, ymin);
3932 
3933  if( c > 0.0 )
3934  {
3935  SCIP_Real sqrtc;
3936 
3937  sqrtc = sqrt(c);
3938  minvalleft = MIN(-sqrtc - b, minvalleft);
3939  maxvalright = MAX( sqrtc - b, maxvalright);
3940  }
3941  }
3942  }
3943  }
3944 
3945  if( rhs.inf > -infinity )
3946  {
3947  sqrtterm = axy * axy * ay * (ay * bx * bx - axy * bx * by + ax * by * by - axy * axy * rhs.inf + 4.0 * ax * ay * rhs.inf);
3948  if( !EPSN(sqrtterm, 1e-9) )
3949  {
3950  sqrtterm = sqrt(MAX(sqrtterm, 0.0));
3951  /* check first candidate for extreme points of +/-sqrt(r(rhs,y))-b(y) */
3952  ymin = axy * ay * bx - 2.0 * ax * ay * by - sqrtterm;
3953  ymin /= ay;
3954  ymin /= 4.0 * ax * ay - axy * axy;
3955 
3956  if( ymin > ybnds.inf && ymin < ybnds.sup )
3957  {
3958  SCIP_Real b;
3959  SCIP_Real c;
3960 
3961  b = CALCB(ymin);
3962  c = CALCR(rhs.inf, ymin);
3963 
3964  if( c > 0.0 )
3965  {
3966  SCIP_Real sqrtc;
3967 
3968  sqrtc = sqrt(c);
3969  maxvalleft = MAX(-sqrtc - b, maxvalleft);
3970  minvalright = MIN( sqrtc - b, minvalright);
3971  }
3972  }
3973 
3974  /* check second candidate for extreme points of +/-sqrt(c(y))-b(y) */
3975  ymin = axy * ay * bx - 2.0 * ax * ay * by + sqrtterm;
3976  ymin /= ay;
3977  ymin /= 4.0 * ax * ay - axy * axy;
3978 
3979  if( ymin > ybnds.inf && ymin < ybnds.sup )
3980  {
3981  SCIP_Real b;
3982  SCIP_Real c;
3983 
3984  b = CALCB(ymin);
3985  c = CALCR(rhs.inf, ymin);
3986 
3987  if( c > 0.0 )
3988  {
3989  SCIP_Real sqrtc;
3990 
3991  sqrtc = sqrt(c);
3992  maxvalleft = MAX(-sqrtc - b, maxvalleft);
3993  minvalright = MIN( sqrtc - b, minvalright);
3994  }
3995  }
3996  }
3997  }
3998  }
3999  else if( REALABS(2.0 * ay * bx - axy * by) > 1e-9 )
4000  {
4001  if( rhs.sup < infinity )
4002  {
4003  ymin = - (4.0 * ay * bx * by - axy * by * by + 4.0 * axy * ay * rhs.sup);
4004  ymin /= 4.0 * ay;
4005  ymin /= 2.0 * ay * bx - axy * by;
4006 
4007  if( ymin > ybnds.inf && ymin < ybnds.sup )
4008  {
4009  SCIP_Real b;
4010  SCIP_Real c;
4011 
4012  b = CALCB(ymin);
4013  c = CALCR(rhs.sup, ymin);
4014 
4015  if( c > 0.0 )
4016  {
4017  SCIP_Real sqrtc;
4018 
4019  sqrtc = sqrt(c);
4020  minvalleft = MIN(-sqrtc - b, minvalleft);
4021  maxvalright = MAX( sqrtc - b, maxvalright);
4022  }
4023  }
4024  }
4025 
4026  if( rhs.inf > -infinity )
4027  {
4028  ymin = - (4.0 * ay * bx * by - axy * by * by + 4.0 * axy * ay * rhs.inf);
4029  ymin /= 4.0 * ay;
4030  ymin /= 2.0 * ay * bx - axy * by;
4031 
4032  if( ymin > ybnds.inf && ymin < ybnds.sup )
4033  {
4034  SCIP_Real b;
4035  SCIP_Real c;
4036 
4037  b = CALCB(ymin);
4038  c = CALCR(rhs.inf, ymin);
4039 
4040  if( c > 0.0 )
4041  {
4042  SCIP_Real sqrtc;
4043 
4044  sqrtc = sqrt(c);
4045  maxvalleft = MAX(-sqrtc - b, maxvalleft);
4046  minvalright = MIN( sqrtc - b, minvalright);
4047  }
4048  }
4049  }
4050  }
4051  }
4052 
4053  /* evaluate the case r(rhs,y) = 0, which is to min/max -b(y) w.r.t. r(rhs,y) = 0, y in ybnds
4054  * with the above assignments
4055  * rcoef_y = axy * bx / (2.0*ax) - by;
4056  * rcoef_yy = axy * axy / (4.0*ax) - ay;
4057  * rcoef_const = bx * bx / (4.0*ax);
4058  * we have r(rhs,y) = rhs + rcoef_const + rcoef_y * y + rcoef_yy * y^2
4059  *
4060  * thus, r(rhs,y) = 0 <-> rcoef_y * y + rcoef_yy * y^2 in -rhs - rcoef_const
4061  *
4062  */
4063  {
4064  SCIP_INTERVAL rcoef_yy_int;
4065  SCIP_INTERVAL rcoef_y_int;
4066  SCIP_INTERVAL rhs2;
4067  SCIP_Real b;
4068 
4069  /* setup rcoef_yy, rcoef_y and -rhs-rcoef_const as intervals */
4070  SCIPintervalSet(&rcoef_yy_int, (SCIP_Real)rcoef_yy);
4071  SCIPintervalSet(&rcoef_y_int, (SCIP_Real)rcoef_y);
4072  SCIPintervalSetBounds(&rhs2, (SCIP_Real)(-rhs.sup - rcoef_const), (SCIP_Real)(-rhs.inf - rcoef_const));
4073 
4074  /* first find all y >= 0 such that rcoef_y * y + rcoef_yy * y^2 in -rhs2, if ybnds.sup >= 0.0
4075  * and evaluate -b(y) w.r.t. these values
4076  */
4077  if( ybnds.sup >= 0.0 )
4078  {
4079  SCIP_INTERVAL ypos;
4080 
4081  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, &ypos, rcoef_yy_int, rcoef_y_int, rhs2, ybnds);
4082  if( !SCIPintervalIsEmpty(infinity, ypos) )
4083  {
4084  assert(ypos.inf >= 0.0); /* we computed only positive solutions above */
4085  b = CALCB(ypos.inf);
4086  minvalleft = MIN(minvalleft, -b);
4087  maxvalleft = MAX(maxvalleft, -b);
4088  minvalright = MIN(minvalright, -b);
4089  maxvalright = MAX(maxvalright, -b);
4090 
4091  if( ypos.sup < infinity )
4092  {
4093  b = CALCB(ypos.sup);
4094  minvalleft = MIN(minvalleft, -b);
4095  maxvalleft = MAX(maxvalleft, -b);
4096  minvalright = MIN(minvalright, -b);
4097  maxvalright = MAX(maxvalright, -b);
4098  }
4099  else
4100  {
4101  /* -b(y) = - (bx + axy * y) / (2.0 * sqrt(ax)) -> -sign(axy)*infinity for y -> infinity */
4102  if( axy > 0.0 )
4103  {
4104  minvalleft = -infinity;
4105  minvalright = -infinity;
4106  }
4107  else
4108  {
4109  maxvalleft = infinity;
4110  maxvalright = infinity;
4111  }
4112  }
4113  }
4114  }
4115 
4116  /* next find all y <= 0 such that rcoef_y * y + rcoef_yy * y^2 in -rhs2, if ybnds.inf < 0.0
4117  * and evaluate -b(y) w.r.t. these values
4118  * (the case y fixed to 0 has been handled in the ybnds.sup >= 0 case above)
4119  */
4120  if( ybnds.inf < 0.0 )
4121  {
4122  SCIP_INTERVAL yneg;
4123 
4124  SCIPintervalSolveUnivariateQuadExpressionNegative(infinity, &yneg, rcoef_yy_int, rcoef_y_int, rhs2, ybnds);
4125  if( !SCIPintervalIsEmpty(infinity, yneg) )
4126  {
4127  if( yneg.inf > -infinity )
4128  {
4129  b = CALCB(yneg.inf);
4130  minvalleft = MIN(minvalleft, -b);
4131  maxvalleft = MAX(maxvalleft, -b);
4132  minvalright = MIN(minvalright, -b);
4133  maxvalright = MAX(maxvalright, -b);
4134  }
4135  else
4136  {
4137  /* -b(y) = - (bx + axy * y) / (2.0 * sqrt(ax)) -> sign(axy)*infinity for y -> -infinity */
4138  if( axy > 0.0 )
4139  {
4140  maxvalleft = infinity;
4141  maxvalright = infinity;
4142  }
4143  else
4144  {
4145  minvalleft = -infinity;
4146  minvalright = -infinity;
4147  }
4148  }
4149 
4150  assert(yneg.sup <= 0.0); /* we computed only negative solutions above */
4151  b = CALCB(yneg.sup);
4152  minvalleft = MIN(minvalleft, -b);
4153  maxvalleft = MAX(maxvalleft, -b);
4154  minvalright = MIN(minvalright, -b);
4155  maxvalright = MAX(maxvalright, -b);
4156  }
4157  }
4158  }
4159 
4160  if( rhs.inf > -infinity && xbnds.inf > -infinity && EPSGT(xbnds.inf, maxvalleft / sqrtax, 1e-9) )
4161  {
4162  /* if sqrt(ax)*x > -sqrt(r(rhs,y))-b(y), then tighten lower bound of sqrt(ax)*x to lower bound of sqrt(r(rhs,y))-b(y)
4163  * this is only possible if rhs.inf > -infinity, otherwise the value for maxvalleft is not valid (but tightening wouldn't be possible for sure anyway) */
4164  assert(EPSGE(minvalright, minvalleft, 1e-9)); /* right interval should not be above lower bound of left interval */
4165  if( minvalright > -infinity )
4166  {
4167  assert(minvalright < infinity);
4168  resultant->inf = (SCIP_Real)(minvalright / sqrtax);
4169  }
4170  }
4171  else
4172  {
4173  /* otherwise, tighten lower bound of sqrt(ax)*x to lower bound of -sqrt(r(rhs,y))-b(y) */
4174  if( minvalleft > -infinity )
4175  {
4176  assert(minvalleft < infinity);
4177  resultant->inf = (SCIP_Real)(minvalleft / sqrtax);
4178  }
4179  }
4180 
4181  if( rhs.inf > -infinity && xbnds.sup < infinity && EPSLT(xbnds.sup, minvalright / sqrtax, 1e-9) )
4182  {
4183  /* if sqrt(ax)*x < sqrt(r(rhs,y))-b(y), then tighten upper bound of sqrt(ax)*x to upper bound of -sqrt(r(rhs,y))-b(y)
4184  * this is only possible if rhs.inf > -infinity, otherwise the value for minvalright is not valid (but tightening wouldn't be possible for sure anyway) */
4185  assert(EPSLE(maxvalleft, maxvalright, 1e-9)); /* left interval should not be above upper bound of right interval */
4186  if( maxvalleft < infinity )
4187  {
4188  assert(maxvalleft > -infinity);
4189  resultant->sup = (SCIP_Real)(maxvalleft / sqrtax);
4190  }
4191  }
4192  else
4193  {
4194  /* otherwise, tighten upper bound of sqrt(ax)*x to upper bound of sqrt(r(rhs,y))-b(y) */
4195  if( maxvalright < infinity )
4196  {
4197  assert(maxvalright > -infinity);
4198  resultant->sup = (SCIP_Real)(maxvalright / sqrtax);
4199  }
4200  }
4201 
4202  resultant->inf -= 1e-10 * REALABS(resultant->inf);
4203  resultant->sup += 1e-10 * REALABS(resultant->sup);
4204 
4205 #undef CALCB
4206 #undef CALCR
4207  }
4208  else
4209  {
4210  /* case ax == 0 */
4211 
4212  SCIP_Real c;
4213  SCIP_Real d;
4214  SCIP_Real ymin;
4215  SCIP_Real minval;
4216  SCIP_Real maxval;
4217 
4218  /* consider -bx / axy in ybnds, i.e., bx + axy y can be 0 */
4219  if( EPSGE(-bx / axy, ybnds.inf, 1e-9) && EPSLE(-bx / axy, ybnds.sup, 1e-9) )
4220  {
4221  /* write as (bx + axy y) * x \in (c - ay y^2 - by y)
4222  * and estimate bx + axy y and c - ay y^2 - by y by intervals independently
4223  * @todo can we do better, as in the case where bx + axy y is bounded away from 0?
4224  */
4225  SCIP_INTERVAL lincoef;
4226  SCIP_INTERVAL myrhs;
4227  SCIP_INTERVAL tmp;
4228 
4229  if( xbnds.inf < 0.0 && xbnds.sup > 0.0 )
4230  {
4231  /* if (bx + axy y) can be arbitrary small and x be both positive and negative,
4232  * then nothing we can tighten here
4233  */
4234  SCIPintervalSetBounds(resultant, xbnds.inf, xbnds.sup);
4235  return;
4236  }
4237 
4238  /* store interval for (bx + axy y) in lincoef */
4239  SCIPintervalMulScalar(infinity, &lincoef, ybnds, axy);
4240  SCIPintervalAddScalar(infinity, &lincoef, lincoef, bx);
4241 
4242  /* store interval for (c - ay y^2 - by y) in myrhs */
4243  SCIPintervalSet(&tmp, by);
4244  SCIPintervalQuad(infinity, &tmp, ay, tmp, ybnds);
4245  SCIPintervalSub(infinity, &myrhs, rhs, tmp);
4246 
4247  if( lincoef.inf == 0.0 && lincoef.sup == 0.0 )
4248  {
4249  /* equation became 0.0 \in myrhs */
4250  if( myrhs.inf <= 0.0 && myrhs.sup >= 0.0 ) /*lint !e644*/
4251  *resultant = xbnds;
4252  else
4253  SCIPintervalSetEmpty(resultant);
4254  }
4255  else if( xbnds.inf >= 0.0 )
4256  {
4257  SCIP_INTERVAL a_;
4258 
4259  /* need only positive solutions */
4260  SCIPintervalSet(&a_, 0.0);
4261  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, resultant, a_, lincoef, myrhs, xbnds);
4262  }
4263  else
4264  {
4265  SCIP_INTERVAL a_;
4266  SCIP_INTERVAL xbndsneg;
4267 
4268  assert(xbnds.sup <= 0.0);
4269 
4270  /* need only negative solutions */
4271  SCIPintervalSet(&a_, 0.0);
4272  SCIPintervalSetBounds(&lincoef, -lincoef.sup, -lincoef.inf);
4273  SCIPintervalSetBounds(&xbndsneg, -xbnds.sup, -xbnds.inf);
4274  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, resultant, a_, lincoef, myrhs, xbndsneg);
4275  if( !SCIPintervalIsEmpty(infinity, *resultant) )
4276  SCIPintervalSetBounds(resultant, -resultant->sup, -resultant->inf);
4277  }
4278 
4279  return;
4280  }
4281 
4282  minval = infinity;
4283  maxval = -infinity;
4284 
4285  /* compute a lower bound on x */
4286  if( bx + axy * (axy > 0.0 ? ybnds.inf : ybnds.sup) > 0.0 )
4287  c = rhs.inf;
4288  else
4289  c = rhs.sup;
4290 
4291  if( c > -infinity && c < infinity )
4292  {
4293  if( ybnds.inf <= -infinity )
4294  {
4295  /* limit is ay/axy * infinity if ay != 0.0 and -by/axy otherwise */
4296  if( EPSZ(ay, 1e-9) )
4297  minval = -by / axy;
4298  else if( ay * axy < 0.0 )
4299  minval = -infinity;
4300  }
4301  else
4302  {
4303  val = (c - ay * ybnds.inf * ybnds.inf - by * ybnds.inf) / (bx + axy * ybnds.inf);
4304  minval = MIN(val, minval);
4305  }
4306 
4307  if( ybnds.sup >= infinity )
4308  {
4309  /* limit is -ay/axy * infinity if ay != 0.0 and -by/axy otherwise */
4310  if( EPSZ(ay, 1e-9) )
4311  minval = MIN(minval, -by / axy);
4312  else if( ay * axy > 0.0 )
4313  minval = -infinity;
4314  }
4315  else
4316  {
4317  val = (c - ay * ybnds.sup * ybnds.sup - by * ybnds.sup) / (bx + axy * ybnds.sup);
4318  minval = MIN(val, minval);
4319  }
4320 
4321  if( !EPSZ(ay, 1e-9) )
4322  {
4323  d = ay * (ay * bx * bx - axy * (bx * by + axy * c));
4324  if( !EPSN(d, 1e-9) )
4325  {
4326  ymin = -ay * bx + sqrt(MAX(d, 0.0));
4327  ymin /= axy * ay;
4328 
4329  if( ymin > ybnds.inf && ymin < ybnds.sup )
4330  {
4331  assert(bx + axy * ymin != 0.0);
4332 
4333  val = (c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4334  minval = MIN(val, minval);
4335  }
4336 
4337  ymin = -ay * bx - sqrt(MAX(d, 0.0));
4338  ymin /= axy * ay;
4339 
4340  if(ymin > ybnds.inf && ymin < ybnds.sup )
4341  {
4342  assert(bx + axy * ymin != 0.0);
4343 
4344  val = (c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4345  minval = MIN(val, minval);
4346  }
4347  }
4348  }
4349  }
4350  else
4351  {
4352  minval = -infinity;
4353  }
4354 
4355  /* compute an upper bound on x */
4356  if( bx + axy * (axy > 0.0 ? ybnds.inf : ybnds.sup) > 0.0 )
4357  c = rhs.sup;
4358  else
4359  c = rhs.inf;
4360 
4361  if( c > -infinity && c < infinity )
4362  {
4363  if( ybnds.inf <= -infinity )
4364  {
4365  /* limit is ay/axy * infinity if ay != 0.0 and -by/axy otherwise */
4366  if( EPSZ(ay, 1e-9) )
4367  maxval = -by / axy;
4368  else if( ay * axy > 0.0 )
4369  maxval = infinity;
4370  }
4371  else
4372  {
4373  val = (c - ay * ybnds.inf * ybnds.inf - by * ybnds.inf) / (bx + axy * ybnds.inf);
4374  maxval = MAX(val, maxval);
4375  }
4376 
4377  if( ybnds.sup >= infinity )
4378  {
4379  /* limit is -ay/axy * infinity if ay != 0.0 and -by/axy otherwise */
4380  if( EPSZ(ay, 1e-9) )
4381  maxval = MAX(maxval, -by / axy);
4382  else if( ay * axy < 0.0 )
4383  maxval = infinity;
4384  }
4385  else
4386  {
4387  val = (c - ay * ybnds.sup * ybnds.sup - by * ybnds.sup) / (bx + axy * ybnds.sup);
4388  maxval = MAX(val, maxval);
4389  }
4390 
4391  if( !EPSZ(ay, 1e-9) )
4392  {
4393  d = ay * (ay * bx * bx - axy * (bx * by + axy * c));
4394  if( !EPSN(d, 1e-9) )
4395  {
4396  ymin = ay * bx + sqrt(MAX(d, 0.0));
4397  ymin /= axy * ay;
4398 
4399  if( ymin > ybnds.inf && ymin < ybnds.sup )
4400  {
4401  assert(bx + axy * ymin != 0.0); /* the case -bx/axy in ybnds was handled aboved */
4402  val = (c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4403  maxval = MAX(val, maxval);
4404  }
4405 
4406  ymin = ay * bx - sqrt(MAX(d, 0.0));
4407  ymin /= axy * ay;
4408 
4409  if( ymin > ybnds.inf && ymin < ybnds.sup )
4410  {
4411  assert(bx + axy * ymin != 0.0); /* the case -bx/axy in ybnds was handled aboved */
4412  val = (c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4413  maxval = MAX(val, maxval);
4414  }
4415  }
4416  }
4417  }
4418  else
4419  {
4420  maxval = infinity;
4421  }
4422 
4423  if( minval > -infinity )
4424  resultant->inf = minval - 1e-10 * REALABS(minval);
4425  else
4426  resultant->inf = -infinity;
4427  if( maxval < infinity )
4428  resultant->sup = maxval + 1e-10 * REALABS(maxval);
4429  else
4430  resultant->sup = infinity;
4431  SCIPintervalIntersect(resultant, *resultant, xbnds);
4432  }
4433 }
4434 
4435 /* pop -O0 from beginning, though it probably doesn't matter here at the end of the compilation unit */
4436 #if defined(__GNUC__) && !defined( __INTEL_COMPILER)
4437 #pragma GCC pop_options
4438 #endif
void SCIPintervalSignPowerScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalDivScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalMulSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSubScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalMax(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalSign(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalAddSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
#define SCIP_ROUND_NEAREST
void SCIPintervalSolveUnivariateQuadExpressionNegative(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
SCIP_Real SCIPintervalNegateReal(SCIP_Real x)
#define infinity
Definition: gastrans.c:71
SCIPInterval pow(const SCIPInterval &x, const SCIPInterval &y)
void SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real sqrcoeff, SCIP_Real lincoeff, SCIP_Real rhs, SCIP_INTERVAL xbnds)
SCIPInterval cos(const SCIPInterval &x)
void SCIPintervalSetRoundingMode(SCIP_ROUNDMODE roundmode)
#define FALSE
Definition: def.h:73
void SCIPintervalMul(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
#define EPSISINT(x, eps)
Definition: def.h:200
#define TRUE
Definition: def.h:72
SCIP_Real SCIPintervalPowerScalarIntegerInf(SCIP_Real operand1, int operand2)
void SCIPintervalMin(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIPInterval exp(const SCIPInterval &x)
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
void SCIPintervalPowerScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
static SCIP_Real negate(SCIP_Real x)
#define EPSGE(x, y, eps)
Definition: def.h:192
#define SCIPdebugMessage
Definition: pub_message.h:87
SCIP_Bool SCIPintervalIsNegativeInfinity(SCIP_Real infinity, SCIP_INTERVAL operand)
SCIP_Real SCIPnegateReal(SCIP_Real x)
Definition: misc.c:10194
SCIP_VAR ** x
Definition: circlepacking.c:54
void SCIPintervalDiv(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
static SCIP_ROUNDMODE intervalGetRoundingMode(void)
void SCIPintervalAddVectors(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_INTERVAL *operand2)
SCIP_Bool SCIPintervalIsPositiveInfinity(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalSin(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
SCIP_Real inf
Definition: intervalarith.h:39
#define EPSN(x, eps)
Definition: def.h:195
void SCIPintervalPowerScalarScalar(SCIP_INTERVAL *resultant, SCIP_Real operand1, SCIP_Real operand2)
SCIP_Bool SCIPintervalIsEntire(SCIP_Real infinity, SCIP_INTERVAL operand)
#define SCIPerrorMessage
Definition: pub_message.h:55
void SCIPintervalScalprod(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_INTERVAL *operand2)
interval arithmetics for provable bounds
#define SCIPdebugPrintf
Definition: pub_message.h:90
void SCIPintervalSetEmpty(SCIP_INTERVAL *resultant)
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
void SCIPintervalLog(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
SCIPInterval sqrt(const SCIPInterval &x)
void SCIPintervalScalprodScalarsInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_Real *operand2)
static SCIP_Real sin_extremepoints[]
Definition: intervalarith.c:42
#define SCIP_ROUND_ZERO
internal miscellaneous methods
#define NULL
Definition: lpi_spx1.cpp:155
#define REALABS(x)
Definition: def.h:187
#define SCIP_ROUND_UPWARDS
SCIP_Real sup
Definition: intervalarith.h:40
void SCIPintervalSet(SCIP_INTERVAL *resultant, SCIP_Real value)
void SCIPintervalPowerScalarInteger(SCIP_INTERVAL *resultant, SCIP_Real operand1, int operand2)
void SCIPintervalCos(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
#define SCIP_Bool
Definition: def.h:70
void SCIPintervalSolveUnivariateQuadExpression(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
SCIPInterval sin(const SCIPInterval &x)
static SCIP_Real cos_extremepoints[]
Definition: intervalarith.c:44
void SCIPintervalSolveBivariateQuadExpressionAllScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real ax, SCIP_Real ay, SCIP_Real axy, SCIP_Real bx, SCIP_Real by, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds, SCIP_INTERVAL ybnds)
#define MAX(x, y)
Definition: tclique_def.h:83
void SCIPintervalSquareRoot(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalQuadBivar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real ax, SCIP_Real ay, SCIP_Real axy, SCIP_Real bx, SCIP_Real by, SCIP_INTERVAL xbnds, SCIP_INTERVAL ybnds)
void SCIPintervalSquare(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalScalprodScalars(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_Real *operand2)
SCIP_Bool SCIPintervalHasRoundingControl(void)
void SCIPintervalSetRoundingModeTowardsZero(void)
#define EPSLE(x, y, eps)
Definition: def.h:190
void SCIPintervalAdd(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIP_Real SCIPintervalGetSup(SCIP_INTERVAL interval)
SCIP_Bool SCIPintervalAreDisjoint(SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIPInterval log(const SCIPInterval &x)
void SCIPintervalMulScalarSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
#define SCIP_ROUND_DOWNWARDS
void SCIPintervalAbs(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalMulScalarInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_Real SCIPnextafter(SCIP_Real from, SCIP_Real to)
Definition: misc.c:9265
#define M_PI
Definition: pricer_rpa.c:88
#define SCIP_REAL_MAX
Definition: def.h:164
#define EPSLT(x, y, eps)
Definition: def.h:189
#define SCIP_REAL_MIN
Definition: def.h:165
#define EPSGT(x, y, eps)
Definition: def.h:191
SCIP_VAR ** b
Definition: circlepacking.c:56
void SCIPintervalSetEntire(SCIP_Real infinity, SCIP_INTERVAL *resultant)
SCIP_Real SCIPintervalQuadUpperBound(SCIP_Real infinity, SCIP_Real a, SCIP_INTERVAL b_, SCIP_INTERVAL x)
void SCIPintervalSetRoundingModeUpwards(void)
void SCIPintervalExp(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
public methods for message output
#define CALCB(y)
SCIP_VAR * a
Definition: circlepacking.c:57
SCIP_Real SCIPintervalPowerScalarIntegerSup(SCIP_Real operand1, int operand2)
#define SCIP_Real
Definition: def.h:163
void SCIPintervalIntersect(SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalMulScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_VAR ** y
Definition: circlepacking.c:55
void SCIPintervalSolveUnivariateQuadExpressionPositive(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
static void intervalSetRoundingMode(SCIP_ROUNDMODE roundmode)
void SCIPintervalSetRoundingModeToNearest(void)
int SCIP_ROUNDMODE
Definition: intervalarith.h:46
void SCIPintervalMulInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSub(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
#define CALCR(c, y)
common defines and data types used in all packages of SCIP
void SCIPintervalAddScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_ROUNDMODE SCIPintervalGetRoundingMode(void)
void SCIPintervalPower(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSetRoundingModeDownwards(void)
#define EPSZ(x, eps)
Definition: def.h:193
void SCIPintervalReciprocal(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalPowerScalarInverse(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL basedomain, SCIP_Real exponent, SCIP_INTERVAL image)
SCIP_Bool SCIPintervalIsSubsetEQ(SCIP_Real infinity, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalAddInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalScalprodScalarsSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_Real *operand2)
void SCIPintervalUnify(SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalQuad(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL xrng)