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