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