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-2023 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24 
25 /**@file intervalarith.c
26  * @ingroup OTHER_CFILES
27  * @brief interval arithmetics for provable bounds
28  * @author Tobias Achterberg
29  * @author Stefan Vigerske
30  * @author Kati Wolter
31  */
32 
33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
34 
35 #define _USE_MATH_DEFINES /* to get M_PI on Windows */ /*lint !750 */
36 
37 #include <stdlib.h>
38 #include <assert.h>
39 #include <math.h>
40 
41 #include "scip/def.h"
42 #include "scip/intervalarith.h"
43 #include "scip/pub_message.h"
44 #include "scip/misc.h"
45 
46 /* Inform compiler that this code accesses the floating-point environment, so that
47  * certain optimizations should be omitted (http://www.cplusplus.com/reference/cfenv/FENV_ACCESS/).
48  * Not supported by Clang (gives warning) and GCC (silently), at the moment.
49  */
50 #if defined(__INTEL_COMPILER) || defined(_MSC_VER)
51 #pragma fenv_access (on)
52 #elif defined(__GNUC__) && !defined(__clang__)
53 #pragma STDC FENV_ACCESS ON
54 #endif
55 
56 /* Unfortunately, the FENV_ACCESS pragma is essentially ignored by GCC at the moment (2019),
57  * see #2650 and https://gcc.gnu.org/bugzilla/show_bug.cgi?id=34678.
58  * There are ways to work around this by declaring variables volatile or inserting more assembler code,
59  * but there is always the danger that something would be overlooked.
60  * A more drastic but safer way seems to be to just disable all compiler optimizations for this file.
61  * The Intel compiler seems to implement FENV_ACCESS correctly, but also defines __GNUC__.
62  */
63 #if defined(__GNUC__) && !defined(__INTEL_COMPILER)
64 #if defined(__clang__)
65 #pragma clang optimize off
66 #else
67 #pragma GCC optimize ("O0")
68 #endif
69 #endif
70 
71 /*lint -e644*/
72 /*lint -e777*/
73 
74 #ifdef SCIP_ROUNDING_FE
75 #define ROUNDING
76 /*
77  * Linux rounding operations
78  */
79 
80 #include <fenv.h>
81 
82 /** Linux rounding mode settings */
83 #define SCIP_ROUND_DOWNWARDS FE_DOWNWARD /**< round always down */
84 #define SCIP_ROUND_UPWARDS FE_UPWARD /**< round always up */
85 #define SCIP_ROUND_NEAREST FE_TONEAREST /**< round always to nearest */
86 #define SCIP_ROUND_ZERO FE_TOWARDZERO /**< round always towards 0.0 */
87 
88 /** returns whether rounding mode control is available */
90  void
91  )
92 {
93  return TRUE;
94 }
95 
96 /** sets rounding mode of floating point operations */
97 static
99  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
100  )
101 {
102 #ifndef NDEBUG
103  if( fesetround(roundmode) != 0 )
104  {
105  SCIPerrorMessage("error setting rounding mode to %d\n", roundmode);
106  abort();
107  }
108 #else
109  (void) fesetround(roundmode);
110 #endif
111 }
112 
113 /** gets current rounding mode of floating point operations */
114 static
116  void
117  )
118 {
119  return (SCIP_ROUNDMODE)fegetround();
120 }
121 
122 #endif
123 
124 
125 
126 #ifdef SCIP_ROUNDING_FP
127 #define ROUNDING
128 /*
129  * OSF rounding operations
130  */
131 
132 #include <float.h>
133 
134 /** OSF rounding mode settings */
135 #define SCIP_ROUND_DOWNWARDS FP_RND_RM /**< round always down */
136 #define SCIP_ROUND_UPWARDS FP_RND_RP /**< round always up */
137 #define SCIP_ROUND_NEAREST FP_RND_RN /**< round always to nearest */
138 #define SCIP_ROUND_ZERO FP_RND_RZ /**< round always towards 0.0 */
139 
140 /** returns whether rounding mode control is available */
142  void
143  )
144 {
145  return TRUE;
146 }
147 
148 /** sets rounding mode of floating point operations */
149 static
151  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
152  )
153 {
154 #ifndef NDEBUG
155  if( write_rnd(roundmode) != 0 )
156  {
157  SCIPerrorMessage("error setting rounding mode to %d\n", roundmode);
158  abort();
159  }
160 #else
161  (void) write_rnd(roundmode);
162 #endif
163 }
164 
165 /** gets current rounding mode of floating point operations */
166 static
168  void
169  )
170 {
171  return read_rnd();
172 }
173 
174 #endif
175 
176 
177 
178 #ifdef SCIP_ROUNDING_MS
179 #define ROUNDING
180 /*
181  * Microsoft compiler rounding operations
182  */
183 
184 #include <float.h>
185 
186 /** Microsoft rounding mode settings */
187 #define SCIP_ROUND_DOWNWARDS RC_DOWN /**< round always down */
188 #define SCIP_ROUND_UPWARDS RC_UP /**< round always up */
189 #define SCIP_ROUND_NEAREST RC_NEAR /**< round always to nearest */
190 #define SCIP_ROUND_ZERO RC_CHOP /**< round always towards zero */
191 
192 /** returns whether rounding mode control is available */
194  void
195  )
196 {
197  return TRUE;
198 }
199 
200 /** sets rounding mode of floating point operations */
201 static
203  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
204  )
205 {
206 #ifndef NDEBUG
207  if( (_controlfp(roundmode, _MCW_RC) & _MCW_RC) != roundmode )
208  {
209  SCIPerrorMessage("error setting rounding mode to %x\n", roundmode);
210  abort();
211  }
212 #else
213  (void) _controlfp(roundmode, _MCW_RC);
214 #endif
215 }
216 
217 /** gets current rounding mode of floating point operations */
218 static
220  void
221  )
222 {
223  return _controlfp(0, 0) & _MCW_RC;
224 }
225 #endif
226 
227 
228 
229 #ifndef ROUNDING
230 /*
231  * rouding operations not available
232  */
233 #define SCIP_ROUND_DOWNWARDS 0 /**< round always down */
234 #define SCIP_ROUND_UPWARDS 1 /**< round always up */
235 #define SCIP_ROUND_NEAREST 2 /**< round always to nearest */
236 #define SCIP_ROUND_ZERO 3 /**< round always towards zero */
237 
238 /** returns whether rounding mode control is available */
240  void
241  )
242 {
243  return FALSE;
244 }
245 
246 /** sets rounding mode of floating point operations */ /*lint -e715*/
247 static
249  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
250  )
251 { /*lint --e{715}*/
252  SCIPerrorMessage("setting rounding mode not available - interval arithmetic is invalid!\n");
253 }
254 
255 /** gets current rounding mode of floating point operations */
256 static
258  void
259  )
260 {
261  return SCIP_ROUND_NEAREST;
262 }
263 #else
264 #undef ROUNDING
265 #endif
266 
267 /** sets rounding mode of floating point operations */
269  SCIP_ROUNDMODE roundmode /**< rounding mode to activate */
270  )
271 {
272  intervalSetRoundingMode(roundmode);
273 }
274 
275 /** gets current rounding mode of floating point operations */
277  void
278  )
279 {
280  return intervalGetRoundingMode();
281 }
282 
283 #if defined(__GNUC__) && (defined(__i386__) || defined(__x86_64__)) /* gcc or icc compiler on x86 32bit or 64bit */
284 
285 /** gets the negation of a double
286  * Do this in a way that the compiler does not "optimize" it away, which usually does not considers rounding modes.
287  * However, compiling with -frounding-math would allow to return -x here.
288  * @todo We now set the FENV_ACCESS pragma to on, which is the same as -frounding-math, so we might be able to eliminate this.
289  */
290 static
291 double negate(
292  /* we explicitly use double here, since I'm not sure the assembler code would work as it for other float's */
293  double x /**< number that should be negated */
294  )
295 {
296  /* The following line of code is taken from GAOL, http://sourceforge.net/projects/gaol. */
297  __asm volatile ("fldl %1; fchs; fstpl %0" : "=m" (x) : "m" (x));
298  return x;
299 }
300 
301 /* cl or icl compiler on 32bit windows or icl compiler on 64bit windows
302  * cl on 64bit windows does not seem to support inline assembler
303  */
304 #elif defined(_MSC_VER) && (defined(__INTEL_COMPILER) || !defined(_M_X64))
305 
306 /** gets the negation of a double
307  * Do this in a way that the compiler does not "optimize" it away, which usually does not considers rounding modes.
308  */
309 static
310 double negate(
311  /* we explicitly use double here, since I'm not sure the assembler code would work as it for other float's */
312  double x /**< number that should be negated */
313  )
314 {
315  /* The following lines of code are taken from GAOL, http://sourceforge.net/projects/gaol. */
316  __asm {
317  fld x
318  fchs
319  fstp x
320  }
321  return x;
322 }
323 
324 #else /* unknown compiler or MSVS 64bit */
325 
326 /** gets the negation of a double
327  *
328  * Fallback implementation that calls the negation method from misc.o.
329  * Having the implementation in a different object file will hopefully prevent
330  * it from being "optimized away".
331  */
332 static
334  SCIP_Real x /**< number that should be negated */
335  )
336 {
337  return SCIPnegateReal(x);
338 }
339 
340 #endif
341 
342 
343 /** sets rounding mode of floating point operations to downwards rounding */
345  void
346  )
347 {
349 }
350 
351 /** sets rounding mode of floating point operations to upwards rounding */
353  void
354  )
355 {
357 }
358 
359 /** sets rounding mode of floating point operations to nearest rounding */
361  void
362  )
363 {
365 }
366 
367 /** sets rounding mode of floating point operations to towards zero rounding */
369  void
370  )
371 {
373 }
374 
375 /** negates a number in a way that the compiler does not optimize it away */
377  SCIP_Real x /**< number to negate */
378  )
379 {
380  return negate((double)x);
381 }
382 
383 /*
384  * Interval arithmetic operations
385  */
386 
387 /* In debug mode, the following methods are implemented as function calls to ensure
388  * type validity.
389  * In optimized mode, the methods are implemented as defines to improve performance.
390  * However, we want to have them in the library anyways, so we have to undef the defines.
391  */
392 
393 #undef SCIPintervalGetInf
394 #undef SCIPintervalGetSup
395 #undef SCIPintervalSet
396 #undef SCIPintervalSetBounds
397 #undef SCIPintervalSetEmpty
398 #undef SCIPintervalIsEmpty
399 #undef SCIPintervalSetEntire
400 #undef SCIPintervalIsEntire
401 #undef SCIPintervalIsPositiveInfinity
402 #undef SCIPintervalIsNegativeInfinity
403 
404 /** returns infimum of interval */
406  SCIP_INTERVAL interval /**< interval */
407  )
408 {
409  return interval.inf;
410 }
411 
412 /** returns supremum of interval */
414  SCIP_INTERVAL interval /**< interval */
415  )
416 {
417  return interval.sup;
418 }
419 
420 /** stores given value as interval */
422  SCIP_INTERVAL* resultant, /**< interval to store value into */
423  SCIP_Real value /**< value to store */
424  )
425 {
426  assert(resultant != NULL);
427 
428  resultant->inf = value;
429  resultant->sup = value;
430 }
431 
432 /** stores given infimum and supremum as interval */
434  SCIP_INTERVAL* resultant, /**< interval to store value into */
435  SCIP_Real inf, /**< value to store as infimum */
436  SCIP_Real sup /**< value to store as supremum */
437  )
438 {
439  assert(resultant != NULL);
440  assert(inf <= sup);
441 
442  resultant->inf = inf;
443  resultant->sup = sup;
444 }
445 
446 /** sets interval to empty interval, which will be [1.0, -1.0] */
448  SCIP_INTERVAL* resultant /**< resultant interval of operation */
449  )
450 {
451  assert(resultant != NULL);
452 
453  resultant->inf = 1.0;
454  resultant->sup = -1.0;
455 }
456 
457 /** indicates whether interval is empty, i.e., whether inf > sup */
459  SCIP_Real infinity, /**< value for infinity */
460  SCIP_INTERVAL operand /**< operand of operation */
461  )
462 {
463  if( operand.sup >= infinity || operand.inf <= -infinity )
464  return FALSE;
465 
466  return operand.sup < operand.inf;
467 }
468 
469 /** sets interval to entire [-infinity, +infinity] */
471  SCIP_Real infinity, /**< value for infinity */
472  SCIP_INTERVAL* resultant /**< resultant interval of operation */
473  )
474 {
475  assert(resultant != NULL);
476 
477  resultant->inf = -infinity;
478  resultant->sup = infinity;
479 }
480 
481 /** indicates whether interval is entire, i.e., whether inf &le; -infinity and sup &ge; infinity */
483  SCIP_Real infinity, /**< value for infinity */
484  SCIP_INTERVAL operand /**< operand of operation */
485  )
486 {
487  return operand.inf <= -infinity && operand.sup >= infinity;
488 }
489 
490 /** indicates whether interval is positive infinity, i.e., [infinity, infinity] */
492  SCIP_Real infinity, /**< value for infinity */
493  SCIP_INTERVAL operand /**< operand of operation */
494  )
495 {
496  return operand.inf >= infinity && operand.sup >= operand.inf;
497 }
498 
499 /** indicates whether interval is negative infinity, i.e., [-infinity, -infinity] */
501  SCIP_Real infinity, /**< value for infinity */
502  SCIP_INTERVAL operand /**< operand of operation */
503  )
504 {
505  return operand.sup <= -infinity && operand.inf <= operand.sup;
506 }
507 
508 /** indicates whether operand1 is contained in operand2 */
510  SCIP_Real infinity, /**< value for infinity */
511  SCIP_INTERVAL operand1, /**< first operand of operation */
512  SCIP_INTERVAL operand2 /**< second operand of operation */
513  )
514 {
515  /* the empty interval is contained everywhere */
516  if( operand1.inf > operand1.sup )
517  return TRUE;
518 
519  /* something not-empty is not contained in the empty interval */
520  if( operand2.inf > operand2.sup )
521  return FALSE;
522 
523  return (MAX(-infinity, operand1.inf) >= operand2.inf) &&
524  ( MIN( infinity, operand1.sup) <= operand2.sup);
525 }
526 
527 /** indicates whether operand1 and operand2 are disjoint */
529  SCIP_INTERVAL operand1, /**< first operand of operation */
530  SCIP_INTERVAL operand2 /**< second operand of operation */
531  )
532 {
533  return (operand1.sup < operand2.inf) || (operand2.sup < operand1.inf);
534 }
535 
536 /** indicates whether operand1 and operand2 are disjoint with epsilon tolerance
537  *
538  * Returns whether minimal (relative) distance of intervals is larger than epsilon.
539  * Same as `SCIPintervalIsEmpty(SCIPintervalIntersectEps(operand1, operand2))`.
540  */
542  SCIP_Real eps, /**< epsilon */
543  SCIP_INTERVAL operand1, /**< first operand of operation */
544  SCIP_INTERVAL operand2 /**< second operand of operation */
545  )
546 {
547  if( operand1.sup < operand2.inf )
548  return SCIPrelDiff(operand2.inf, operand1.sup) > eps;
549 
550  if( operand1.inf > operand2.sup )
551  return SCIPrelDiff(operand1.inf, operand2.sup) > eps;
552 
553  return FALSE;
554 }
555 
556 /** intersection of two intervals */
558  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
559  SCIP_INTERVAL operand1, /**< first operand of operation */
560  SCIP_INTERVAL operand2 /**< second operand of operation */
561  )
562 {
563  assert(resultant != NULL);
564 
565  resultant->inf = MAX(operand1.inf, operand2.inf);
566  resultant->sup = MIN(operand1.sup, operand2.sup);
567 }
568 
569 /** intersection of two intervals with epsilon tolerance
570  *
571  * If intersection of operand1 and operand2 is empty, but minimal (relative) distance of intervals
572  * is at most epsilon, then set resultant to singleton containing the point in operand1
573  * that is closest to operand2, i.e.,
574  * - `resultant = { operand1.sup }`, if `operand1.sup` < `operand2.inf` and `reldiff(operand2.inf,operand1.sup)` &le; eps
575  * - `resultant = { operand1.inf }`, if `operand1.inf` > `operand2.sup` and `reldiff(operand1.inf,operand2.sup)` &le; eps
576  * - `resultant` = intersection of `operand1` and `operand2`, otherwise
577  */
579  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
580  SCIP_Real eps, /**< epsilon */
581  SCIP_INTERVAL operand1, /**< first operand of operation */
582  SCIP_INTERVAL operand2 /**< second operand of operation */
583  )
584 {
585  assert(resultant != NULL);
586  assert(eps >= 0.0);
587 
588  if( operand1.sup < operand2.inf )
589  {
590  if( SCIPrelDiff(operand2.inf, operand1.sup) <= eps )
591  {
592  SCIPintervalSet(resultant, operand1.sup);
593  return;
594  }
595  }
596  else if( operand1.inf > operand2.sup )
597  {
598  if( SCIPrelDiff(operand1.inf, operand2.sup) <= eps )
599  {
600  SCIPintervalSet(resultant, operand1.inf);
601  return;
602  }
603  }
604 
605  SCIPintervalIntersect(resultant, operand1, operand2);
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 = intervalGetRoundingMode();
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  intervalSetRoundingMode(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 = intervalGetRoundingMode();
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  intervalSetRoundingMode(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 = intervalGetRoundingMode();
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  intervalSetRoundingMode(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 = intervalGetRoundingMode();
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 */
833  resultant->sup = -infinity;
834  }
835  else
836  {
838  resultant->sup = operand1.sup - operand2.inf;
839  }
840 
841  intervalSetRoundingMode(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 = intervalGetRoundingMode();
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  intervalSetRoundingMode(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  if( operand2 == 1.0 )
1128  {
1129  *resultant = operand1;
1130  return;
1131  }
1132 
1133  if( operand2 == -1.0 )
1134  {
1135  resultant->inf = -operand1.sup;
1136  resultant->sup = -operand1.inf;
1137  return;
1138  }
1139 
1140  roundmode = intervalGetRoundingMode();
1141 
1142  /* compute infimum result */
1144  SCIPintervalMulScalarInf(infinity, resultant, operand1, operand2);
1145 
1146  /* compute supremum of result */
1148  SCIPintervalMulScalarSup(infinity, resultant, operand1, operand2);
1149 
1150  intervalSetRoundingMode(roundmode);
1151 }
1152 
1153 /** divides operand1 by operand2 and stores result in resultant */
1155  SCIP_Real infinity, /**< value for infinity */
1156  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1157  SCIP_INTERVAL operand1, /**< first operand of operation */
1158  SCIP_INTERVAL operand2 /**< second operand of operation */
1159  )
1160 {
1161  SCIP_ROUNDMODE roundmode;
1162  SCIP_INTERVAL intmed;
1163 
1164  assert(resultant != NULL);
1165  assert(!SCIPintervalIsEmpty(infinity, operand1));
1166  assert(!SCIPintervalIsEmpty(infinity, operand2));
1167 
1168  if( operand2.inf <= 0.0 && operand2.sup >= 0.0 )
1169  { /* division by [0,0] or interval containing 0 gives [-inf, +inf] */
1170  resultant->inf = -infinity;
1171  resultant->sup = infinity;
1172  return;
1173  }
1174 
1175  if( operand1.inf == 0.0 && operand1.sup == 0.0 )
1176  { /* division of [0,0] by something nonzero */
1177  SCIPintervalSet(resultant, 0.0);
1178  return;
1179  }
1180 
1181  roundmode = intervalGetRoundingMode();
1182 
1183  /* division by nonzero: resultant = x * (1/y) */
1184  if( operand2.sup >= infinity || operand2.sup <= -infinity )
1185  {
1186  intmed.inf = 0.0;
1187  }
1188  else
1189  {
1191  intmed.inf = 1.0 / operand2.sup;
1192  }
1193  if( operand2.inf <= -infinity || operand2.inf >= infinity )
1194  {
1195  intmed.sup = 0.0;
1196  }
1197  else
1198  {
1200  intmed.sup = 1.0 / operand2.inf;
1201  }
1202  SCIPintervalMul(infinity, resultant, operand1, intmed);
1203 
1204  intervalSetRoundingMode(roundmode);
1205 }
1206 
1207 /** divides operand1 by scalar operand2 and stores result in resultant */
1209  SCIP_Real infinity, /**< value for infinity */
1210  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1211  SCIP_INTERVAL operand1, /**< first operand of operation */
1212  SCIP_Real operand2 /**< second operand of operation */
1213  )
1214 {
1215  SCIP_ROUNDMODE roundmode;
1216 
1217  assert(resultant != NULL);
1218  assert(!SCIPintervalIsEmpty(infinity, operand1));
1219 
1220  roundmode = intervalGetRoundingMode();
1221 
1222  if( operand2 >= infinity || operand2 <= -infinity )
1223  {
1224  /* division by +/-infinity is 0.0 */
1225  resultant->inf = 0.0;
1226  resultant->sup = 0.0;
1227  }
1228  else if( operand2 > 0.0 )
1229  {
1230  if( operand1.inf <= -infinity )
1231  resultant->inf = -infinity;
1232  else if( operand1.inf >= infinity )
1233  {
1234  /* infinity / + = infinity */
1235  resultant->inf = infinity;
1236  }
1237  else
1238  {
1240  resultant->inf = operand1.inf / operand2;
1241  }
1242 
1243  if( operand1.sup >= infinity )
1244  resultant->sup = infinity;
1245  else if( operand1.sup <= -infinity )
1246  {
1247  /* -infinity / + = -infinity */
1248  resultant->sup = -infinity;
1249  }
1250  else
1251  {
1253  resultant->sup = operand1.sup / operand2;
1254  }
1255  }
1256  else if( operand2 < 0.0 )
1257  {
1258  if( operand1.sup >= infinity )
1259  resultant->inf = -infinity;
1260  else if( operand1.sup <= -infinity )
1261  {
1262  /* -infinity / - = infinity */
1263  resultant->inf = infinity;
1264  }
1265  else
1266  {
1268  resultant->inf = operand1.sup / operand2;
1269  }
1270 
1271  if( operand1.inf <= -infinity )
1272  resultant->sup = infinity;
1273  else if( operand1.inf >= infinity )
1274  {
1275  /* infinity / - = -infinity */
1276  resultant->sup = -infinity;
1277  }
1278  else
1279  {
1281  resultant->sup = operand1.inf / operand2;
1282  }
1283  }
1284  else
1285  { /* division by 0.0 */
1286  if( operand1.inf >= 0 )
1287  {
1288  /* [+,+] / [0,0] = [+inf, +inf] */
1289  resultant->inf = infinity;
1290  resultant->sup = infinity;
1291  }
1292  else if( operand1.sup <= 0 )
1293  {
1294  /* [-,-] / [0,0] = [-inf, -inf] */
1295  resultant->inf = -infinity;
1296  resultant->sup = -infinity;
1297  }
1298  else
1299  {
1300  /* [-,+] / [0,0] = [-inf, +inf] */
1301  resultant->inf = -infinity;
1302  resultant->sup = infinity;
1303  }
1304  return;
1305  }
1306 
1307  intervalSetRoundingMode(roundmode);
1308 }
1309 
1310 /** computes the scalar product of two vectors of intervals and stores result in resultant */
1312  SCIP_Real infinity, /**< value for infinity */
1313  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1314  int length, /**< length of vectors */
1315  SCIP_INTERVAL* operand1, /**< first vector as array of intervals; can have +/-inf entries */
1316  SCIP_INTERVAL* operand2 /**< second vector as array of intervals; can have +/-inf entries */
1317  )
1318 {
1319  SCIP_ROUNDMODE roundmode;
1320  SCIP_INTERVAL prod;
1321  int i;
1322 
1323  roundmode = intervalGetRoundingMode();
1324 
1325  resultant->inf = 0.0;
1326  resultant->sup = 0.0;
1327 
1328  /* compute infimum of resultant */
1330  for( i = 0; i < length && resultant->inf > -infinity; ++i )
1331  {
1332  SCIPintervalSetEntire(infinity, &prod);
1333  SCIPintervalMulInf(infinity, &prod, operand1[i], operand2[i]);
1334  SCIPintervalAddInf(infinity, resultant, *resultant, prod);
1335  }
1336  assert(resultant->sup == 0.0);
1337 
1338  /* compute supremum of resultant */
1340  for( i = 0; i < length && resultant->sup < infinity ; ++i )
1341  {
1342  SCIPintervalSetEntire(infinity, &prod);
1343  SCIPintervalMulSup(infinity, &prod, operand1[i], operand2[i]);
1344  SCIPintervalAddSup(infinity, resultant, *resultant, prod);
1345  }
1346 
1347  intervalSetRoundingMode(roundmode);
1348 }
1349 
1350 /** computes the scalar product of a vector of intervals and a vector of scalars and stores infimum of result in infimum of resultant */
1352  SCIP_Real infinity, /**< value for infinity */
1353  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1354  int length, /**< length of vectors */
1355  SCIP_INTERVAL* operand1, /**< first vector as array of intervals */
1356  SCIP_Real* operand2 /**< second vector as array of scalars; can have +/-inf entries */
1357  )
1358 {
1359  SCIP_INTERVAL prod;
1360  int i;
1361 
1363 
1364  resultant->inf = 0.0;
1365 
1366  /* compute infimum of resultant */
1367  SCIPintervalSetEntire(infinity, &prod);
1368  for( i = 0; i < length && resultant->inf > -infinity; ++i )
1369  {
1370  SCIPintervalMulScalarInf(infinity, &prod, operand1[i], operand2[i]);
1371  assert(prod.sup >= infinity);
1372  SCIPintervalAddInf(infinity, resultant, *resultant, prod);
1373  }
1374 }
1375 
1376 /** computes the scalar product of a vector of intervals and a vector of scalars and stores supremum of result in supremum of resultant */
1378  SCIP_Real infinity, /**< value for infinity */
1379  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1380  int length, /**< length of vectors */
1381  SCIP_INTERVAL* operand1, /**< first vector as array of intervals */
1382  SCIP_Real* operand2 /**< second vector as array of scalars; can have +/-inf entries */
1383  )
1384 {
1385  SCIP_INTERVAL prod;
1386  int i;
1387 
1389 
1390  resultant->sup = 0.0;
1391 
1392  /* compute supremum of resultant */
1393  SCIPintervalSetEntire(infinity, &prod);
1394  for( i = 0; i < length && resultant->sup < infinity; ++i )
1395  {
1396  SCIPintervalMulScalarSup(infinity, &prod, operand1[i], operand2[i]);
1397  assert(prod.inf <= -infinity);
1398  SCIPintervalAddSup(infinity, resultant, *resultant, prod);
1399  }
1400 }
1401 
1402 /** computes the scalar product of a vector of intervals and a vector of scalars and stores result in resultant */
1404  SCIP_Real infinity, /**< value for infinity */
1405  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1406  int length, /**< length of vectors */
1407  SCIP_INTERVAL* operand1, /**< first vector as array of intervals */
1408  SCIP_Real* operand2 /**< second vector as array of scalars; can have +/-inf entries */
1409  )
1410 {
1411  SCIP_ROUNDMODE roundmode;
1412 
1413  roundmode = intervalGetRoundingMode();
1414 
1415  resultant->inf = 0.0;
1416  resultant->sup = 0.0;
1417 
1418  /* compute infimum of resultant */
1420  SCIPintervalScalprodScalarsInf(infinity, resultant, length, operand1, operand2);
1421  assert(resultant->sup == 0.0);
1422 
1423  /* compute supremum of resultant */
1425  SCIPintervalScalprodScalarsSup(infinity, resultant, length, operand1, operand2);
1426 
1427  intervalSetRoundingMode(roundmode);
1428 }
1429 
1430 /** squares operand and stores result in resultant */
1432  SCIP_Real infinity, /**< value for infinity */
1433  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1434  SCIP_INTERVAL operand /**< operand of operation */
1435  )
1436 {
1437  SCIP_ROUNDMODE roundmode;
1438 
1439  assert(resultant != NULL);
1440  assert(!SCIPintervalIsEmpty(infinity, operand));
1441 
1442  roundmode = intervalGetRoundingMode();
1443 
1444  if( operand.sup <= 0.0 )
1445  { /* operand is left of 0.0 */
1446  if( operand.sup <= -infinity )
1447  resultant->inf = infinity;
1448  else
1449  {
1451  resultant->inf = operand.sup * operand.sup;
1452  }
1453 
1454  if( operand.inf <= -infinity )
1455  resultant->sup = infinity;
1456  else
1457  {
1459  resultant->sup = operand.inf * operand.inf;
1460  }
1461  }
1462  else if( operand.inf >= 0.0 )
1463  { /* operand is right of 0.0 */
1464  if( operand.inf >= infinity )
1465  resultant->inf = infinity;
1466  else
1467  {
1469  resultant->inf = operand.inf * operand.inf;
1470  }
1471 
1472  if( operand.sup >= infinity )
1473  resultant->sup = infinity;
1474  else
1475  {
1477  resultant->sup = operand.sup * operand.sup;
1478  }
1479  }
1480  else
1481  { /* [-,+]^2 */
1482  resultant->inf = 0.0;
1483  if( operand.inf <= -infinity || operand.sup >= infinity )
1484  resultant->sup = infinity;
1485  else
1486  {
1487  SCIP_Real x;
1488  SCIP_Real y;
1489 
1491  x = operand.inf * operand.inf;
1492  y = operand.sup * operand.sup;
1493  resultant->sup = MAX(x, y);
1494  }
1495  }
1496 
1497  intervalSetRoundingMode(roundmode);
1498 }
1499 
1500 /** stores (positive part of) square root of operand in resultant
1501  * @attention we assume a correctly rounded sqrt(double) function when rounding is to nearest
1502  */
1504  SCIP_Real infinity, /**< value for infinity */
1505  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1506  SCIP_INTERVAL operand /**< operand of operation */
1507  )
1508 {
1509  assert(resultant != NULL);
1510  assert(!SCIPintervalIsEmpty(infinity, operand));
1511 
1512  if( operand.sup < 0.0 )
1513  {
1514  SCIPintervalSetEmpty(resultant);
1515  return;
1516  }
1517 
1518  if( operand.inf == operand.sup )
1519  {
1520  if( operand.inf >= infinity )
1521  {
1522  resultant->inf = infinity;
1523  resultant->sup = infinity;
1524  }
1525  else
1526  {
1527  SCIP_Real tmp;
1528 
1529  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1530  tmp = sqrt(operand.inf);
1531  resultant->inf = SCIPnextafter(tmp, SCIP_REAL_MIN);
1532  resultant->sup = SCIPnextafter(tmp, SCIP_REAL_MAX);
1533  }
1534 
1535  return;
1536  }
1537 
1538  if( operand.inf <= 0.0 )
1539  resultant->inf = 0.0;
1540  else if( operand.inf >= infinity )
1541  {
1542  resultant->inf = infinity;
1543  resultant->sup = infinity;
1544  }
1545  else
1546  {
1547  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1548  resultant->inf = SCIPnextafter(sqrt(operand.inf), SCIP_REAL_MIN);
1549  }
1550 
1551  if( operand.sup >= infinity )
1552  resultant->sup = infinity;
1553  else
1554  {
1555  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1556  resultant->sup = SCIPnextafter(sqrt(operand.sup), SCIP_REAL_MAX);
1557  }
1558 }
1559 
1560 /** stores operand1 to the power of operand2 in resultant
1561  *
1562  * uses SCIPintervalPowerScalar if operand2 is a scalar, otherwise computes exp(op2*log(op1))
1563  */
1565  SCIP_Real infinity, /**< value for infinity */
1566  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1567  SCIP_INTERVAL operand1, /**< first operand of operation */
1568  SCIP_INTERVAL operand2 /**< second operand of operation */
1569  )
1570 {
1571  assert(resultant != NULL);
1572  assert(!SCIPintervalIsEmpty(infinity, operand1));
1573  assert(!SCIPintervalIsEmpty(infinity, operand2));
1574 
1575  if( operand2.inf == operand2.sup )
1576  { /* operand is number */
1577  SCIPintervalPowerScalar(infinity, resultant, operand1, operand2.inf);
1578  return;
1579  }
1580 
1581  /* log([..,0]) will give an empty interval below, but we want [0,0]^exponent to be 0
1582  * if 0 is in exponent, then resultant should also contain 1 (the case exponent == [0,0] is handled above)
1583  */
1584  if( operand1.sup == 0.0 )
1585  {
1586  if( operand2.inf <= 0.0 && operand2.sup >= 0.0 )
1587  SCIPintervalSetBounds(resultant, 0.0, 1.0);
1588  else
1589  SCIPintervalSet(resultant, 0.0);
1590  return;
1591  }
1592 
1593  /* resultant := log(op1) */
1594  SCIPintervalLog(infinity, resultant, operand1);
1595  if( SCIPintervalIsEmpty(infinity, *resultant) )
1596  return;
1597 
1598  /* resultant := op2 * resultant */
1599  SCIPintervalMul(infinity, resultant, operand2, *resultant);
1600 
1601  /* resultant := exp(resultant) */
1602  SCIPintervalExp(infinity, resultant, *resultant);
1603 }
1604 
1605 /** computes lower bound on power of a scalar operand1 to an integer operand2
1606  *
1607  * Both operands need to be finite numbers.
1608  * Needs to have operand1 &ge; 0 and need to have operand2 &ge; 0 if operand1 = 0.
1609  */
1611  SCIP_Real operand1, /**< first operand of operation */
1612  int operand2 /**< second operand of operation */
1613  )
1614 {
1615  SCIP_ROUNDMODE roundmode;
1616  SCIP_Real result;
1617 
1618  assert(operand1 >= 0.0);
1619 
1620  if( operand1 == 0.0 )
1621  {
1622  assert(operand2 >= 0);
1623  if( operand2 == 0 )
1624  return 1.0; /* 0^0 = 1 */
1625  else
1626  return 0.0; /* 0^positive = 0 */
1627  }
1628 
1629  /* 1^n = 1, x^0 = 1 */
1630  if( operand1 == 1.0 || operand2 == 0 )
1631  return 1.0;
1632 
1633  if( operand2 < 0 )
1634  {
1635  /* x^n = 1 / x^(-n) */
1636  result = SCIPintervalPowerScalarIntegerSup(operand1, -operand2);
1637  assert(result != 0.0);
1638 
1639  roundmode = intervalGetRoundingMode();
1641  result = 1.0 / result;
1642  intervalSetRoundingMode(roundmode);
1643  }
1644  else
1645  {
1646  unsigned int n;
1647  SCIP_Real z;
1648 
1649  roundmode = intervalGetRoundingMode();
1650 
1651  result = 1.0;
1652  n = (unsigned int)operand2;
1653  z = operand1;
1654 
1656 
1657  /* use a binary exponentiation algorithm:
1658  * consider the binary representation of n: n = sum_i 2^i d_i with d_i \in {0,1}
1659  * then x^n = prod_{i:d_i=1} x^(2^i)
1660  * in the following, we loop over i=1,..., thereby storing x^(2^i) in z
1661  * whenever d_i is 1, we multiply result with x^(2^i) (the current z)
1662  * at the last (highest) i with d_i = 1 we stop, thus having x^n stored in result
1663  *
1664  * the binary representation of n and bit shifting is used for the loop
1665  */
1666  assert(n >= 1);
1667  do
1668  {
1669  if( n & 1 ) /* n is odd (d_i=1), so multiply result with current z (=x^{2^i}) */
1670  {
1671  result = result * z;
1672  n >>= 1;
1673  if( n == 0 )
1674  break;
1675  }
1676  else
1677  n >>= 1;
1678  z = z * z;
1679  }
1680  while( TRUE ); /*lint !e506 */
1681 
1682  intervalSetRoundingMode(roundmode);
1683  }
1684 
1685  return result;
1686 }
1687 
1688 /** computes upper bound on power of a scalar operand1 to an integer operand2
1689  *
1690  * Both operands need to be finite numbers.
1691  * Needs to have operand1 &ge; 0 and needs to have operand2 &ge; 0 if operand1 = 0.
1692  */
1694  SCIP_Real operand1, /**< first operand of operation */
1695  int operand2 /**< second operand of operation */
1696  )
1697 {
1698  SCIP_ROUNDMODE roundmode;
1699  SCIP_Real result;
1700 
1701  assert(operand1 >= 0.0);
1702 
1703  if( operand1 == 0.0 )
1704  {
1705  assert(operand2 >= 0);
1706  if( operand2 == 0 )
1707  return 1.0; /* 0^0 = 1 */
1708  else
1709  return 0.0; /* 0^positive = 0 */
1710  }
1711 
1712  /* 1^n = 1, x^0 = 1 */
1713  if( operand1 == 1.0 || operand2 == 0 )
1714  return 1.0;
1715 
1716  if( operand2 < 0 )
1717  {
1718  /* x^n = 1 / x^(-n) */
1719  result = SCIPintervalPowerScalarIntegerInf(operand1, -operand2);
1720  assert(result != 0.0);
1721 
1722  roundmode = intervalGetRoundingMode();
1724  result = 1.0 / result;
1725  intervalSetRoundingMode(roundmode);
1726  }
1727  else
1728  {
1729  unsigned int n;
1730  SCIP_Real z;
1731 
1732  roundmode = intervalGetRoundingMode();
1733 
1734  result = 1.0;
1735  n = (unsigned int)operand2;
1736  z = operand1;
1737 
1739 
1740  /* use a binary exponentiation algorithm... see comments in SCIPintervalPowerScalarIntegerInf */
1741  assert(n >= 1);
1742  do
1743  {
1744  if( n&1 )
1745  {
1746  result = result * z;
1747  n >>= 1;
1748  if( n == 0 )
1749  break;
1750  }
1751  else
1752  n >>= 1;
1753  z = z * z;
1754  }
1755  while( TRUE ); /*lint !e506 */
1756 
1757  intervalSetRoundingMode(roundmode);
1758  }
1759 
1760  return result;
1761 }
1762 
1763 /** computes bounds on power of a scalar operand1 to an integer operand2
1764  *
1765  * Both operands need to be finite numbers.
1766  * Needs to have operand1 &ge; 0 and needs to have operand2 &ge; 0 if operand1 = 0.
1767  */
1769  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1770  SCIP_Real operand1, /**< first operand of operation */
1771  int operand2 /**< second operand of operation */
1772  )
1773 {
1774  SCIP_ROUNDMODE roundmode;
1775 
1776  assert(operand1 >= 0.0);
1777 
1778  if( operand1 == 0.0 )
1779  {
1780  assert(operand2 >= 0);
1781  if( operand2 == 0 )
1782  {
1783  SCIPintervalSet(resultant, 1.0); /* 0^0 = 1 */
1784  return;
1785  }
1786  else
1787  {
1788  SCIPintervalSet(resultant, 0.0); /* 0^positive = 0 */
1789  return;
1790  }
1791  }
1792 
1793  /* 1^n = 1, x^0 = 1 */
1794  if( operand1 == 1.0 || operand2 == 0 )
1795  {
1796  SCIPintervalSet(resultant, 1.0);
1797  return;
1798  }
1799 
1800  if( operand2 < 0 )
1801  {
1802  /* x^n = 1 / x^(-n) */
1803  SCIPintervalPowerScalarInteger(resultant, operand1, -operand2);
1804  assert(resultant->inf > 0.0 || resultant->sup < 0.0);
1805  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 */
1806  }
1807  else
1808  {
1809  unsigned int n;
1810  SCIP_Real z_inf;
1811  SCIP_Real z_sup;
1812  SCIP_Real result_sup;
1813  SCIP_Real result_inf;
1814 
1815  roundmode = intervalGetRoundingMode();
1816 
1817  result_inf = 1.0;
1818  result_sup = 1.0;
1819  z_inf = operand1;
1820  z_sup = operand1;
1821  n = (unsigned int)operand2;
1822 
1824 
1825  /* use a binary exponentiation algorithm... see comments in SCIPintervalPowerScalarIntegerInf
1826  * we compute lower and upper bounds within the same loop
1827  * to get correct lower bounds while rounding mode is upwards, we negate arguments */
1828  assert(n >= 1);
1829  do
1830  {
1831  if( n & 1 )
1832  {
1833  result_inf = negate(negate(result_inf) * z_inf);
1834  result_sup = result_sup * z_sup;
1835  n >>= 1;
1836  if( n == 0 )
1837  break;
1838  }
1839  else
1840  n >>= 1;
1841  z_inf = negate(negate(z_inf) * z_inf);
1842  z_sup = z_sup * z_sup;
1843  }
1844  while( TRUE ); /*lint !e506 */
1845 
1846  intervalSetRoundingMode(roundmode);
1847 
1848  resultant->inf = result_inf;
1849  resultant->sup = result_sup;
1850  assert(resultant->inf <= resultant->sup);
1851  }
1852 }
1853 
1854 /** stores bounds on the power of a scalar operand1 to a scalar operand2 in resultant
1855  *
1856  * Both operands need to be finite numbers.
1857  * Needs to have operand1 &ge; 0 or operand2 integer and needs to have operand2 &ge; 0 if operand1 = 0.
1858  * @attention we assume a correctly rounded pow(double) function when rounding is to nearest
1859  */
1861  SCIP_INTERVAL* resultant, /**< resultant of operation */
1862  SCIP_Real operand1, /**< first operand of operation */
1863  SCIP_Real operand2 /**< second operand of operation */
1864  )
1865 {
1866  SCIP_Real result;
1867 
1868  assert(resultant != NULL);
1869 
1870  if( operand1 == 0.0 )
1871  {
1872  assert(operand2 >= 0);
1873  if( operand2 == 0 )
1874  {
1875  SCIPintervalSet(resultant, 1.0); /* 0^0 = 1 */
1876  return;
1877  }
1878  else
1879  {
1880  SCIPintervalSet(resultant, 0.0); /* 0^positive = 0 */
1881  return;
1882  }
1883  }
1884 
1885  /* 1^n = 1, x^0 = 1 */
1886  if( operand1 == 1.0 || operand2 == 0 )
1887  {
1888  SCIPintervalSet(resultant, 1.0);
1889  return;
1890  }
1891 
1892  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1893  result = pow(operand1, operand2);
1894 
1895  /* to get safe bounds, get the floating point numbers just below and above result */
1896  resultant->inf = SCIPnextafter(result, SCIP_REAL_MIN);
1897  resultant->sup = SCIPnextafter(result, SCIP_REAL_MAX);
1898 }
1899 
1900 /** stores operand1 to the power of the scalar operand2 in resultant
1901  * @attention we assume a correctly rounded pow(double) function when rounding is to nearest
1902  */
1904  SCIP_Real infinity, /**< value for infinity */
1905  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
1906  SCIP_INTERVAL operand1, /**< first operand of operation */
1907  SCIP_Real operand2 /**< second operand of operation */
1908  )
1909 {
1910  SCIP_Bool op2isint;
1911 
1912  assert(resultant != NULL);
1913  assert(!SCIPintervalIsEmpty(infinity, operand1));
1914 
1915  if( operand2 == infinity )
1916  {
1917  /* 0^infinity = 0
1918  * +^infinity = infinity
1919  * -^infinity = -infinity
1920  */
1921  if( operand1.inf < 0.0 )
1922  resultant->inf = -infinity;
1923  else
1924  resultant->inf = 0.0;
1925  if( operand1.sup > 0.0 )
1926  resultant->sup = infinity;
1927  else
1928  resultant->sup = 0.0;
1929  return;
1930  }
1931 
1932  if( operand2 == 0.0 )
1933  { /* special case, since x^0 = 1 for x != 0, but 0^0 = 0 */
1934  if( operand1.inf == 0.0 && operand1.sup == 0.0 )
1935  {
1936  resultant->inf = 0.0;
1937  resultant->sup = 0.0;
1938  }
1939  else if( operand1.inf <= 0.0 || operand1.sup >= 0.0 )
1940  { /* 0.0 in x gives [0,1] */
1941  resultant->inf = 0.0;
1942  resultant->sup = 1.0;
1943  }
1944  else
1945  { /* 0.0 outside x gives [1,1] */
1946  resultant->inf = 1.0;
1947  resultant->sup = 1.0;
1948  }
1949  return;
1950  }
1951 
1952  if( operand2 == 1.0 )
1953  {
1954  /* x^1 = x */
1955  *resultant = operand1;
1956  return;
1957  }
1958 
1959  op2isint = (ceil(operand2) == operand2);
1960 
1961  if( !op2isint && operand1.inf < 0.0 )
1962  { /* x^n with x negative not defined for n not integer*/
1963  operand1.inf = 0.0;
1964  if( operand1.sup < operand1.inf )
1965  {
1966  SCIPintervalSetEmpty(resultant);
1967  return;
1968  }
1969  }
1970 
1971  if( operand1.inf >= 0.0 )
1972  { /* easy case: x^n with x >= 0 */
1973  if( operand2 >= 0.0 )
1974  {
1975  /* inf^+ = inf */
1976  if( operand1.inf >= infinity )
1977  resultant->inf = infinity;
1978  else if( operand1.inf > 0.0 )
1979  {
1980  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1981  resultant->inf = SCIPnextafter(pow(operand1.inf, operand2), SCIP_REAL_MIN);
1982  }
1983  else
1984  resultant->inf = 0.0;
1985 
1986  if( operand1.sup >= infinity )
1987  resultant->sup = infinity;
1988  else if( operand1.sup > 0.0 )
1989  {
1990  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
1991  resultant->sup = SCIPnextafter(pow(operand1.sup, operand2), SCIP_REAL_MAX);
1992  }
1993  else
1994  resultant->sup = 0.0;
1995  }
1996  else
1997  {
1998  if( operand1.sup >= infinity )
1999  resultant->inf = 0.0;
2000  else if( operand1.sup == 0.0 )
2001  {
2002  /* x^(negative even) = infinity for x->0 (from both sides),
2003  * but x^(negative odd) = -infinity for x->0 from left side */
2004  if( ceil(operand2/2) == operand2/2 )
2005  resultant->inf = infinity;
2006  else
2007  resultant->inf = -infinity;
2008  }
2009  else
2010  {
2011  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
2012  resultant->inf = SCIPnextafter(pow(operand1.sup, operand2), SCIP_REAL_MIN);
2013  }
2014 
2015  /* 0^(negative) = infinity */
2016  if( operand1.inf == 0.0 )
2017  resultant->sup = infinity;
2018  else
2019  {
2020  assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST); /* usually, no-one should have changed rounding mode */
2021  resultant->sup = SCIPnextafter(pow(operand1.inf, operand2), SCIP_REAL_MAX);
2022  }
2023  }
2024  }
2025  else if( operand1.sup <= 0.0 )
2026  { /* more difficult case: x^n with x < 0; we now know, that n is integer */
2027  assert(op2isint);
2028  if( operand2 >= 0.0 && ceil(operand2/2) == operand2/2 )
2029  {
2030  /* x^n with n>=2 and even -> x^n is monotonically decreasing for x < 0 */
2031  if( operand1.sup == -infinity )
2032  /* (-inf)^n = inf */
2033  resultant->inf = infinity;
2034  else
2035  resultant->inf = SCIPintervalPowerScalarIntegerInf(-operand1.sup, (int)operand2);
2036 
2037  if( operand1.inf <= -infinity )
2038  resultant->sup = infinity;
2039  else
2040  resultant->sup = SCIPintervalPowerScalarIntegerSup(-operand1.inf, (int)operand2);
2041  }
2042  else if( operand2 <= 0.0 && ceil(operand2/2) != operand2/2 )
2043  {
2044  /* x^n with n<=-1 and odd -> x^n = 1/x^(-n) is monotonically decreasing for x<0 */
2045  if( operand1.sup == -infinity )
2046  /* (-inf)^n = 1/(-inf)^(-n) = 1/(-inf) = 0 */
2047  resultant->inf = 0.0;
2048  else if( operand1.sup == 0.0 )
2049  /* x^n -> -infinity for x->0 from left */
2050  resultant->inf = -infinity;
2051  else
2052  resultant->inf = -SCIPintervalPowerScalarIntegerSup(-operand1.sup, (int)operand2);
2053 
2054  if( operand1.inf <= -infinity )
2055  /* (-inf)^n = 1/(-inf)^(-n) = 1/(-inf) = 0 */
2056  resultant->sup = 0.0;
2057  else if( operand1.inf == 0.0 )
2058  /* x^n -> infinity for x->0 from right */
2059  resultant->sup = infinity;
2060  else
2061  resultant->sup = -SCIPintervalPowerScalarIntegerInf(-operand1.inf, (int)operand2);
2062  }
2063  else if( operand2 >= 0.0 )
2064  {
2065  /* x^n with n>0 and odd -> x^n is monotonically increasing for x<0 */
2066  if( operand1.inf <= -infinity )
2067  resultant->inf = -infinity;
2068  else
2069  resultant->inf = -SCIPintervalPowerScalarIntegerSup(-operand1.inf, (int)operand2);
2070 
2071  if( operand1.sup <= -infinity )
2072  resultant->sup = -infinity;
2073  else
2074  resultant->sup = -SCIPintervalPowerScalarIntegerInf(-operand1.sup, (int)operand2);
2075  }
2076  else
2077  {
2078  /* x^n with n<0 and even -> x^n is monotonically increasing for x<0 */
2079  if( operand1.inf <= -infinity )
2080  resultant->inf = 0.0;
2081  else if( operand1.inf == 0.0 )
2082  /* x^n -> infinity for x->0 from both sides */
2083  resultant->inf = infinity;
2084  else
2085  resultant->inf = SCIPintervalPowerScalarIntegerSup(-operand1.inf, (int)operand2);
2086 
2087  if( operand1.sup <= -infinity )
2088  resultant->sup = 0.0;
2089  else if( operand1.sup == 0.0 )
2090  /* x^n -> infinity for x->0 from both sides */
2091  resultant->sup = infinity;
2092  else
2093  resultant->sup = SCIPintervalPowerScalarIntegerSup(-operand1.sup, (int)operand2);
2094  }
2095  assert(resultant->inf <= resultant->sup || resultant->inf >= infinity || resultant->sup <= -infinity);
2096  }
2097  else
2098  { /* similar difficult case: x^n with x in [<0, >0], but n is integer */
2099  assert(op2isint); /* otherwise we had set operand1.inf == 0.0, which was handled in first case */
2100  if( operand2 >= 0.0 && operand2/2 == ceil(operand2/2) )
2101  {
2102  /* n even positive integer */
2103  resultant->inf = 0.0;
2104  if( operand1.inf == -infinity || operand1.sup == infinity )
2105  resultant->sup = infinity;
2106  else
2107  resultant->sup = SCIPintervalPowerScalarIntegerSup(MAX(-operand1.inf, operand1.sup), (int)operand2);
2108  }
2109  else if( operand2 <= 0.0 && ceil(operand2/2) == operand2/2 )
2110  {
2111  /* n even negative integer */
2112  resultant->sup = infinity; /* since 0^n = infinity */
2113  if( operand1.inf == -infinity || operand1.sup == infinity )
2114  resultant->inf = 0.0;
2115  else
2116  resultant->inf = SCIPintervalPowerScalarIntegerInf(MAX(-operand1.inf, operand1.sup), (int)operand2);
2117  }
2118  else if( operand2 >= 0.0 )
2119  {
2120  /* n odd positive integer, so monotonically increasing function */
2121  if( operand1.inf == -infinity )
2122  resultant->inf = -infinity;
2123  else
2124  resultant->inf = -SCIPintervalPowerScalarIntegerSup(-operand1.inf, (int)operand2);
2125  if( operand1.sup == infinity )
2126  resultant->sup = infinity;
2127  else
2128  resultant->sup = SCIPintervalPowerScalarIntegerSup(operand1.sup, (int)operand2);
2129  }
2130  else
2131  {
2132  /* n odd negative integer:
2133  * x^n -> -infinity for x->0 from left
2134  * x^n -> infinity for x->0 from right */
2135  resultant->inf = -infinity;
2136  resultant->sup = infinity;
2137  }
2138  }
2139 
2140  /* if value for infinity is too small, relax intervals so they do not appear empty */
2141  if( resultant->inf > infinity )
2142  resultant->inf = infinity;
2143  if( resultant->sup < -infinity )
2144  resultant->sup = -infinity;
2145 }
2146 
2147 /** given an interval for the image of a power operation, computes an interval for the origin
2148  *
2149  * That is, for \f$y = x^p\f$ with the exponent \f$p\f$ a given scalar and \f$y\f$ = `image` a given interval,
2150  * computes \f$x \subseteq \text{basedomain}\f$ such that \f$y \in x^p\f$ and such that for all \f$z \in \text{basedomain} \setminus x: z^p \not \in y\f$.
2151  */
2153  SCIP_Real infinity, /**< value for infinity */
2154  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2155  SCIP_INTERVAL basedomain, /**< domain of base */
2156  SCIP_Real exponent, /**< exponent */
2157  SCIP_INTERVAL image /**< interval image of power */
2158  )
2159 {
2160  SCIP_INTERVAL tmp;
2161  SCIP_INTERVAL exprecip;
2162 
2163  assert(resultant != NULL);
2164  assert(image.inf <= image.sup);
2165  assert(basedomain.inf <= basedomain.sup);
2166 
2167  if( exponent == 0.0 )
2168  {
2169  /* exponent is 0.0 */
2170  if( image.inf <= 1.0 && image.sup >= 1.0 )
2171  {
2172  /* 1 in image -> resultant = entire */
2173  *resultant = basedomain;
2174  }
2175  else if( image.inf <= 0.0 && image.sup >= 0.0 )
2176  {
2177  /* 0 in image, 1 not in image -> resultant = 0 (provided 0^0 = 0 ???)
2178  * -> resultant = {0} intersected with basedomain */
2179  SCIPintervalSetBounds(resultant, MAX(0.0, basedomain.inf), MIN(0.0, basedomain.sup));
2180  }
2181  else
2182  {
2183  /* 0 and 1 not in image -> resultant = empty */
2184  SCIPintervalSetEmpty(resultant);
2185  }
2186  return;
2187  }
2188 
2189  /* i = b^e
2190  * i >= 0 -> b = i^(1/e) [union -i^(1/e), if e is even]
2191  * i < 0, e odd integer -> b = -(-i)^(1/e)
2192  * i < 0, e even integer or fractional -> empty
2193  */
2194 
2195  SCIPintervalSetBounds(&exprecip, exponent, exponent);
2196  SCIPintervalReciprocal(infinity, &exprecip, exprecip);
2197 
2198  /* invert positive part of image, if any */
2199  if( image.sup >= 0.0 )
2200  {
2201  SCIPintervalSetBounds(&tmp, MAX(image.inf, 0.0), image.sup);
2202  SCIPintervalPower(infinity, resultant, tmp, exprecip);
2203  if( basedomain.inf <= -resultant->inf && EPSISINT(exponent, 0.0) && (int)exponent % 2 == 0 ) /*lint !e835 */
2204  {
2205  if( basedomain.sup < resultant->inf )
2206  SCIPintervalSetBounds(resultant, -resultant->sup, -resultant->inf);
2207  else
2208  SCIPintervalSetBounds(resultant, -resultant->sup, resultant->sup);
2209  }
2210 
2211  SCIPintervalIntersect(resultant, *resultant, basedomain);
2212  }
2213  else
2214  SCIPintervalSetEmpty(resultant);
2215 
2216  /* invert negative part of image, if any and if base can take negative value and if exponent is such that negative values are possible */
2217  if( image.inf < 0.0 && basedomain.inf < 0.0 && EPSISINT(exponent, 0.0) && ((int)exponent % 2 != 0) ) /*lint !e835 */
2218  {
2219  SCIPintervalSetBounds(&tmp, MAX(-image.sup, 0.0), -image.inf);
2220  SCIPintervalPower(infinity, &tmp, tmp, exprecip);
2221  SCIPintervalSetBounds(&tmp, -tmp.sup, -tmp.inf);
2222  SCIPintervalIntersect(&tmp, basedomain, tmp);
2223  SCIPintervalUnify(resultant, *resultant, tmp);
2224  }
2225 }
2226 
2227 /** stores operand1 to the signed power of the scalar positive operand2 in resultant
2228  *
2229  * The signed power of x w.r.t. an exponent n &ge; 0 is given as \f$\mathrm{sign}(x) |x|^n\f$.
2230  *
2231  * @attention we assume correctly rounded sqrt(double) and pow(double) functions when rounding is to nearest
2232  */
2234  SCIP_Real infinity, /**< value for infinity */
2235  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2236  SCIP_INTERVAL operand1, /**< first operand of operation */
2237  SCIP_Real operand2 /**< second operand of operation */
2238  )
2239 {
2240  SCIP_ROUNDMODE roundmode;
2241  assert(resultant != NULL);
2242 
2243  assert(!SCIPintervalIsEmpty(infinity, operand1));
2244  assert(operand2 >= 0.0);
2245 
2246  if( operand2 == infinity )
2247  {
2248  /* 0^infinity = 0
2249  * +^infinity = infinity
2250  *-+^infinity = -infinity
2251  */
2252  if( operand1.inf < 0.0 )
2253  resultant->inf = -infinity;
2254  else
2255  resultant->inf = 0.0;
2256  if( operand1.sup > 0.0 )
2257  resultant->sup = infinity;
2258  else
2259  resultant->sup = 0.0;
2260  return;
2261  }
2262 
2263  if( operand2 == 0.0 )
2264  {
2265  /* special case, since x^0 = 1 for x != 0, but 0^0 = 0 */
2266  if( operand1.inf < 0.0 )
2267  resultant->inf = -1.0;
2268  else if( operand1.inf == 0.0 )
2269  resultant->inf = 0.0;
2270  else
2271  resultant->inf = 1.0;
2272 
2273  if( operand1.sup < 0.0 )
2274  resultant->sup = -1.0;
2275  else if( operand1.sup == 0.0 )
2276  resultant->sup = 0.0;
2277  else
2278  resultant->sup = 1.0;
2279 
2280  return;
2281  }
2282 
2283  if( operand2 == 1.0 )
2284  { /* easy case that should run fast */
2285  *resultant = operand1;
2286  return;
2287  }
2288 
2289  roundmode = intervalGetRoundingMode();
2290 
2291  if( operand2 == 2.0 )
2292  { /* common case where pow can easily be avoided */
2293  if( operand1.inf <= -infinity )
2294  {
2295  resultant->inf = -infinity;
2296  }
2297  else if( operand1.inf >= infinity )
2298  {
2299  resultant->inf = infinity;
2300  }
2301  else if( operand1.inf > 0.0 )
2302  {
2304  resultant->inf = operand1.inf * operand1.inf;
2305  }
2306  else
2307  {
2308  /* need upwards since we negate result of multiplication */
2310  resultant->inf = negate(operand1.inf * operand1.inf);
2311  }
2312 
2313  if( operand1.sup >= infinity )
2314  {
2315  resultant->sup = infinity;
2316  }
2317  else if( operand1.sup <= -infinity )
2318  {
2319  resultant->sup = -infinity;
2320  }
2321  else if( operand1.sup > 0.0 )
2322  {
2324  resultant->sup = operand1.sup * operand1.sup;
2325  }
2326  else
2327  {
2328  /* need downwards since we negate result of multiplication */
2330  resultant->sup = negate(operand1.sup * operand1.sup);
2331  }
2332  assert(resultant->inf <= resultant->sup);
2333  }
2334  else if( operand2 == 0.5 )
2335  { /* another common case where pow can easily be avoided */
2336  if( operand1.inf <= -infinity )
2337  resultant->inf = -infinity;
2338  else if( operand1.inf >= infinity )
2339  resultant->inf = infinity;
2340  else if( operand1.inf >= 0.0 )
2341  {
2343  resultant->inf = SCIPnextafter(sqrt( operand1.inf), SCIP_REAL_MIN);
2344  }
2345  else
2346  {
2348  resultant->inf = -SCIPnextafter(sqrt(-operand1.inf), SCIP_REAL_MAX);
2349  }
2350 
2351  if( operand1.sup >= infinity )
2352  resultant->sup = infinity;
2353  else if( operand1.sup <= -infinity )
2354  resultant->sup = -infinity;
2355  else if( operand1.sup > 0.0 )
2356  {
2358  resultant->sup = SCIPnextafter(sqrt( operand1.sup), SCIP_REAL_MAX);
2359  }
2360  else
2361  {
2363  resultant->sup = -SCIPnextafter(sqrt(-operand1.sup), SCIP_REAL_MAX);
2364  }
2365  assert(resultant->inf <= resultant->sup);
2366  }
2367  else
2368  {
2369  if( operand1.inf <= -infinity )
2370  resultant->inf = -infinity;
2371  else if( operand1.inf >= infinity )
2372  resultant->inf = infinity;
2373  else if( operand1.inf > 0.0 )
2374  {
2376  resultant->inf = SCIPnextafter(pow( operand1.inf, operand2), SCIP_REAL_MIN);
2377  }
2378  else
2379  {
2381  resultant->inf = -SCIPnextafter(pow(-operand1.inf, operand2), SCIP_REAL_MAX);
2382  }
2383 
2384  if( operand1.sup >= infinity )
2385  resultant->sup = infinity;
2386  else if( operand1.sup <= -infinity )
2387  resultant->sup = -infinity;
2388  else if( operand1.sup > 0.0 )
2389  {
2391  resultant->sup = SCIPnextafter(pow( operand1.sup, operand2), SCIP_REAL_MAX);
2392  }
2393  else
2394  {
2396  resultant->sup = -SCIPnextafter(pow(-operand1.sup, operand2), SCIP_REAL_MIN);
2397  }
2398  }
2399 
2400  intervalSetRoundingMode(roundmode);
2401 }
2402 
2403 /** computes the reciprocal of an interval
2404  */
2406  SCIP_Real infinity, /**< value for infinity */
2407  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2408  SCIP_INTERVAL operand /**< operand of operation */
2409  )
2410 {
2411  SCIP_ROUNDMODE roundmode;
2412 
2413  assert(resultant != NULL);
2414  assert(!SCIPintervalIsEmpty(infinity, operand));
2415 
2416  if( operand.inf == 0.0 && operand.sup == 0.0 )
2417  { /* 1/0 = [-inf,inf] */
2418  resultant->inf = infinity;
2419  resultant->sup = -infinity;
2420  return;
2421  }
2422 
2423  roundmode = intervalGetRoundingMode();
2424 
2425  if( operand.inf >= 0.0 )
2426  { /* 1/x with x >= 0 */
2427  if( operand.sup >= infinity )
2428  resultant->inf = 0.0;
2429  else
2430  {
2432  resultant->inf = 1.0 / operand.sup;
2433  }
2434 
2435  if( operand.inf >= infinity )
2436  resultant->sup = 0.0;
2437  else if( operand.inf == 0.0 )
2438  resultant->sup = infinity;
2439  else
2440  {
2442  resultant->sup = 1.0 / operand.inf;
2443  }
2444 
2445  intervalSetRoundingMode(roundmode);
2446  }
2447  else if( operand.sup <= 0.0 )
2448  { /* 1/x with x <= 0 */
2449  if( operand.sup <= -infinity )
2450  resultant->inf = 0.0;
2451  else if( operand.sup == 0.0 )
2452  resultant->inf = -infinity;
2453  else
2454  {
2456  resultant->inf = 1.0 / operand.sup;
2457  }
2458 
2459  if( operand.inf <= -infinity )
2460  resultant->sup = infinity;
2461  else
2462  {
2464  resultant->sup = 1.0 / operand.inf;
2465  }
2466  intervalSetRoundingMode(roundmode);
2467  }
2468  else
2469  { /* 1/x with x in [-,+] is division by zero */
2470  resultant->inf = -infinity;
2471  resultant->sup = infinity;
2472  }
2473 }
2474 
2475 /** stores exponential of operand in resultant
2476  * @attention we assume a correctly rounded exp(double) function when rounding is to nearest
2477  */
2479  SCIP_Real infinity, /**< value for infinity */
2480  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2481  SCIP_INTERVAL operand /**< operand of operation */
2482  )
2483 {
2484  assert(resultant != NULL);
2485  assert(!SCIPintervalIsEmpty(infinity, operand));
2486 
2487  if( operand.sup <= -infinity )
2488  {
2489  resultant->inf = 0.0;
2490  resultant->sup = 0.0;
2491  return;
2492  }
2493 
2494  if( operand.inf >= infinity )
2495  {
2496  resultant->inf = infinity;
2497  resultant->sup = infinity;
2498  return;
2499  }
2500 
2501  if( operand.inf == operand.sup )
2502  {
2503  if( operand.inf == 0.0 )
2504  {
2505  resultant->inf = 1.0;
2506  resultant->sup = 1.0;
2507  }
2508  else
2509  {
2510  SCIP_Real tmp;
2511 
2513  tmp = exp(operand.inf);
2514  resultant->inf = tmp > 0.0 ? SCIPnextafter(tmp, SCIP_REAL_MIN) : 0.0;
2515  assert(resultant->inf >= 0.0);
2516  resultant->sup = SCIPnextafter(tmp, SCIP_REAL_MAX);
2517 
2518  return;
2519  }
2520  }
2521 
2522  if( operand.inf <= -infinity )
2523  {
2524  resultant->inf = 0.0;
2525  }
2526  else if( operand.inf == 0.0 )
2527  {
2528  resultant->inf = 1.0;
2529  }
2530  else
2531  {
2532  SCIP_Real tmp;
2533 
2535  tmp = exp(operand.inf);
2536  resultant->inf = tmp > 0.0 ? SCIPnextafter(tmp, SCIP_REAL_MIN) : 0.0;
2537  /* make sure we do not exceed value for infinity, so interval is not declared as empty if inf and sup are both > infinity */
2538  if( resultant->inf >= infinity )
2539  resultant->inf = infinity;
2540  }
2541 
2542  if( operand.sup >= infinity )
2543  {
2544  resultant->sup = infinity;
2545  }
2546  else if( operand.sup == 0.0 )
2547  {
2548  resultant->sup = 1.0;
2549  }
2550  else
2551  {
2553  resultant->sup = SCIPnextafter(exp(operand.sup), SCIP_REAL_MAX);
2554  if( resultant->sup < -infinity )
2555  resultant->sup = -infinity;
2556  }
2557 }
2558 
2559 /** stores natural logarithm of operand in resultant
2560  * @attention we assume a correctly rounded log(double) function when rounding is to nearest
2561  */
2563  SCIP_Real infinity, /**< value for infinity */
2564  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2565  SCIP_INTERVAL operand /**< operand of operation */
2566  )
2567 {
2568  assert(resultant != NULL);
2569  assert(!SCIPintervalIsEmpty(infinity, operand));
2570 
2571  /* if operand.sup == 0.0, we could return -inf in resultant->sup, but that
2572  * seems of little use and just creates problems somewhere else, e.g., #1230
2573  */
2574  if( operand.sup <= 0.0 )
2575  {
2576  SCIPintervalSetEmpty(resultant);
2577  return;
2578  }
2579 
2580  if( operand.inf == operand.sup )
2581  {
2582  if( operand.sup == 1.0 )
2583  {
2584  resultant->inf = 0.0;
2585  resultant->sup = 0.0;
2586  }
2587  else
2588  {
2589  SCIP_Real tmp;
2590 
2592  tmp = log(operand.inf);
2593  resultant->inf = SCIPnextafter(tmp, SCIP_REAL_MIN);
2594  resultant->sup = SCIPnextafter(tmp, SCIP_REAL_MAX);
2595  }
2596 
2597  return;
2598  }
2599 
2600  if( operand.inf <= 0.0 )
2601  {
2602  resultant->inf = -infinity;
2603  }
2604  else if( operand.inf == 1.0 )
2605  {
2606  resultant->inf = 0.0;
2607  }
2608  else
2609  {
2611  resultant->inf = SCIPnextafter(log(operand.inf), SCIP_REAL_MIN);
2612  }
2613 
2614  if( operand.sup >= infinity )
2615  {
2616  resultant->sup = infinity;
2617  }
2618  else if( operand.sup == 1.0 )
2619  {
2620  resultant->sup = 0.0;
2621  }
2622  else
2623  {
2625  resultant->sup = SCIPnextafter(log(operand.sup), SCIP_REAL_MAX);
2626  }
2627 }
2628 
2629 /** stores minimum of operands in resultant */
2631  SCIP_Real infinity, /**< value for infinity */
2632  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2633  SCIP_INTERVAL operand1, /**< first operand of operation */
2634  SCIP_INTERVAL operand2 /**< second operand of operation */
2635  )
2636 {
2637  assert(resultant != NULL);
2638  assert(!SCIPintervalIsEmpty(infinity, operand1));
2639  assert(!SCIPintervalIsEmpty(infinity, operand2));
2640 
2641  resultant->inf = MIN(operand1.inf, operand2.inf);
2642  resultant->sup = MIN(operand1.sup, operand2.sup);
2643 }
2644 
2645 /** stores maximum of operands in resultant */
2647  SCIP_Real infinity, /**< value for infinity */
2648  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2649  SCIP_INTERVAL operand1, /**< first operand of operation */
2650  SCIP_INTERVAL operand2 /**< second operand of operation */
2651  )
2652 {
2653  assert(resultant != NULL);
2654  assert(!SCIPintervalIsEmpty(infinity, operand1));
2655  assert(!SCIPintervalIsEmpty(infinity, operand2));
2656 
2657  resultant->inf = MAX(operand1.inf, operand2.inf);
2658  resultant->sup = MAX(operand1.sup, operand2.sup);
2659 }
2660 
2661 /** stores absolute value of operand in resultant */
2663  SCIP_Real infinity, /**< value for infinity */
2664  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2665  SCIP_INTERVAL operand /**< operand of operation */
2666  )
2667 {
2668  assert(resultant != NULL);
2669  assert(!SCIPintervalIsEmpty(infinity, operand));
2670 
2671  if( operand.inf <= 0.0 && operand.sup >= 0.0)
2672  {
2673  resultant->inf = 0.0;
2674  resultant->sup = MAX(-operand.inf, operand.sup);
2675  }
2676  else if( operand.inf > 0.0 )
2677  {
2678  *resultant = operand;
2679  }
2680  else
2681  {
2682  resultant->inf = -operand.sup;
2683  resultant->sup = -operand.inf;
2684  }
2685 }
2686 
2687 /* double precision lower and upper bounds on pi
2688  * taken from boost::numeric::interval_lib::constants
2689  * MSVC refuses to evaluate this at compile time
2690  */
2691 #ifndef _MSC_VER
2692 static const double pi_d_l = (3373259426.0 + 273688.0 / (1<<21)) / (1<<30); /*lint !e790*/
2693 static const double pi_d_u = (3373259426.0 + 273689.0 / (1<<21)) / (1<<30); /*lint !e790*/
2694 #else
2695 #define pi_d_l ((3373259426.0 + 273688.0 / (1<<21)) / (1<<30))
2696 #define pi_d_u ((3373259426.0 + 273689.0 / (1<<21)) / (1<<30))
2697 #endif
2698 
2699 /** stores sine value of operand in resultant */
2701  SCIP_Real infinity, /**< value for infinity */
2702  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2703  SCIP_INTERVAL operand /**< operand of operation */
2704  )
2705 {
2706  /* the function evaluates sine transforming it to a cosine via sin(x) = cos(x-pi/2) = -cos(x+pi/2) */
2707  SCIP_INTERVAL pihalf;
2708  SCIP_INTERVAL shiftedop;
2709 
2710  /* sin(x) = cos(x-pi/2) = -cos(x+pi/2)*/
2711  SCIPintervalSetBounds(&pihalf, pi_d_l, pi_d_u);
2712  SCIPintervalMulScalar(infinity, &pihalf, pihalf, 0.5);
2713 
2714  /* intervalCos() will move operand.inf into [0,pi]
2715  * if we can achieve this here by add pi/2 instead of subtracting it, then use the sin(x) = -cos(x+pi/2) identity
2716  */
2717  if( operand.inf < 0.0 && operand.inf > -pi_d_l )
2718  {
2719  SCIP_Real tmp;
2720 
2721  SCIPintervalAdd(infinity, &shiftedop, operand, pihalf);
2722  SCIPintervalCos(infinity, resultant, shiftedop);
2723 
2724  tmp = -resultant->sup;
2725  resultant->sup = -resultant->inf;
2726  resultant->inf = tmp;
2727  }
2728  else
2729  {
2730  SCIPintervalSub(infinity, &shiftedop, operand, pihalf);
2731  SCIPintervalCos(infinity, resultant, shiftedop);
2732  }
2733 
2734  /* some correction if inf or sup is 0, then sin(0) = 0 would be nice */
2735  if( operand.inf == 0.0 && operand.sup < pi_d_l )
2736  resultant->inf = 0.0;
2737  else if( operand.sup == 0.0 && operand.inf > -pi_d_l )
2738  resultant->sup = 0.0;
2739 }
2740 
2741 /** stores cosine value of operand in resultant */
2743  SCIP_Real infinity, /**< value for infinity */
2744  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2745  SCIP_INTERVAL operand /**< operand of operation */
2746  )
2747 {
2748  /* this implementation follows boost::numeric::cos
2749  * cos is decreasing in [0, pi] and increasing in [pi, 2pi].
2750  * If operand = [a,b] and a is in [0, pi], then
2751  * cos([a,b]) = [-1, 1] if b >= 2pi
2752  * cos([a,b]) = [-1, max(cos(a), cos(b))] if b is in [pi, 2pi]
2753  * cos([a,b]) = [cos(b), cos(a)] if b is in [0, pi]
2754  *
2755  * To make sure that a is always between [0, pi] we use the identity cos(x) = (-1)^k cos(x + k pi), i.e.,
2756  * we compute k such that a + k pi \in [0,pi], compute cos([a,b] + k pi) and then multiply by (-1)^k.
2757  */
2758  SCIP_ROUNDMODE roundmode;
2759  SCIP_Real negwidth;
2760  SCIP_Real k = 0.0;
2761 
2762  assert(resultant != NULL);
2763  assert(!SCIPintervalIsEmpty(infinity, operand));
2764 
2765  SCIPdebugMessage("cos([%.16g,%.16g])\n", operand.inf, operand.sup);
2766 
2767  if( operand.inf == operand.sup )
2768  {
2769  SCIP_Real tmp;
2770 
2772  tmp = cos(operand.inf);
2773  resultant->inf = SCIPnextafter(tmp, SCIP_REAL_MIN);
2774  resultant->sup = SCIPnextafter(tmp, SCIP_REAL_MAX);
2775  return;
2776  }
2777 
2778  /* set interval to [-1,1] if we cannot reliably work out the difference between inf and sup
2779  * double precision has almost 16 digits of precision; for now cut off at 12
2780  */
2781  if( operand.sup > 1e12 || operand.inf < -1e12 )
2782  {
2783  SCIPintervalSetBounds(resultant, -1.0, 1.0);
2784  return;
2785  }
2786 
2787  roundmode = SCIPintervalGetRoundingMode();
2788 
2789  /* set interval to [-1,1] if width is at least 2 pi */
2791  negwidth = operand.inf - operand.sup;
2792  if( -negwidth >= 2.0*pi_d_l )
2793  {
2794  SCIPintervalSetBounds(resultant, -1.0, 1.0);
2795  SCIPintervalSetRoundingMode(roundmode);
2796  return;
2797  }
2798 
2799  /* get operand.inf into [0,pi] */
2800  if( operand.inf < 0.0 || operand.inf >= pi_d_l )
2801  {
2802  SCIP_INTERVAL tmp;
2803 
2804  k = floor((operand.inf / (operand.inf < 0.0 ? pi_d_l : pi_d_u)));
2805 
2806  /* operand <- operand - k * pi */
2807  SCIPintervalSetBounds(&tmp, pi_d_l, pi_d_u);
2808  SCIPintervalMulScalar(infinity, &tmp, tmp, k);
2809  SCIPintervalSub(infinity, &operand, operand, tmp);
2810  }
2811  assert(operand.inf >= 0.0);
2812  assert(operand.inf <= pi_d_u);
2813 
2814  SCIPdebugMessage("shifted operand by %g*pi = [%.16g,%.16g])\n", k, operand.inf, operand.sup);
2815 
2816  SCIPintervalSetRoundingMode(roundmode);
2817 
2818  if( operand.sup <= pi_d_l )
2819  {
2820  /* monotone decreasing */
2821  resultant->inf = SCIPnextafter(cos(operand.sup), SCIP_REAL_MIN);
2822  resultant->inf = MAX(-1.0, resultant->inf);
2823  if( operand.inf == 0.0 )
2824  resultant->sup = 1.0;
2825  else
2826  {
2827  resultant->sup = SCIPnextafter(cos(operand.inf), SCIP_REAL_MAX);
2828  resultant->sup = MIN( 1.0, resultant->sup);
2829  }
2830  SCIPdebugMessage("cos([%.16g,%.16g]) = [%.16g,%.16g]\n", operand.inf, operand.sup, resultant->inf, resultant->sup);
2831  }
2832  else if( operand.sup <= 2*pi_d_l )
2833  {
2834  /* inf <= pi, sup >= pi: minimum at pi (=-1), maximum at inf or sup */
2835  resultant->inf = -1.0;
2836  if( operand.inf == 0.0 )
2837  resultant->sup = 1.0;
2838  else
2839  {
2840  SCIP_Real cinf;
2841  SCIP_Real csup;
2842 
2843  cinf = cos(operand.inf);
2844  csup = cos(operand.sup);
2845  resultant->sup = SCIPnextafter(MAX(cinf, csup), SCIP_REAL_MAX);
2846  resultant->sup = MIN(1.0, resultant->sup);
2847  }
2848  SCIPdebugMessage("cos([%.16g,%.16g]) = [%.16g,%.16g]\n", operand.inf, operand.sup, resultant->inf, resultant->sup);
2849  }
2850  else
2851  {
2852  SCIPintervalSetBounds(resultant, -1.0, 1.0);
2853  }
2854 
2855  /* back to original operand using cos(x + k pi) = (-1)^k cos(x) */
2856  if( (int)k % 2 != 0 )
2857  {
2858  SCIP_Real tmp = -resultant->sup;
2859  resultant->sup = -resultant->inf;
2860  resultant->inf = tmp;
2861  SCIPdebugMessage("shifted back -> [%.16g,%.16g]\n", resultant->inf, resultant->sup);
2862  }
2863 
2864  assert(resultant->inf >= -1.0);
2865  assert(resultant->sup <= 1.0);
2866 }
2867 
2868 /** stores sign of operand in resultant */
2870  SCIP_Real infinity, /**< value for infinity */
2871  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2872  SCIP_INTERVAL operand /**< operand of operation */
2873  )
2874 {
2875  assert(resultant != NULL);
2876  assert(!SCIPintervalIsEmpty(infinity, operand));
2877 
2878  if( operand.sup < 0.0 )
2879  {
2880  resultant->inf = -1.0;
2881  resultant->sup = -1.0;
2882  }
2883  else if( operand.inf >= 0.0 )
2884  {
2885  resultant->inf = 1.0;
2886  resultant->sup = 1.0;
2887  }
2888  else
2889  {
2890  resultant->inf = -1.0;
2891  resultant->sup = 1.0;
2892  }
2893 }
2894 
2895 /** stores entropy of operand in resultant */
2897  SCIP_Real infinity, /**< value for infinity */
2898  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
2899  SCIP_INTERVAL operand /**< operand of operation */
2900  )
2901 {
2902  SCIP_Real loginf;
2903  SCIP_Real logsup;
2904  SCIP_Real infcand1 = 0.0;
2905  SCIP_Real infcand2 = 0.0;
2906  SCIP_Real supcand1 = 0.0;
2907  SCIP_Real supcand2 = 0.0;
2908  SCIP_Real extr;
2909  SCIP_Real inf;
2910  SCIP_Real sup;
2911 
2912  assert(resultant != NULL);
2913  assert(!SCIPintervalIsEmpty(infinity, operand));
2914 
2915  /* check whether the domain is empty */
2916  if( operand.sup < 0.0 )
2917  {
2918  SCIPintervalSetEmpty(resultant);
2919  return;
2920  }
2921 
2922  /* handle special case of domain being [0,0] */
2923  if( operand.sup == 0.0 )
2924  {
2925  SCIPintervalSet(resultant, 0.0);
2926  return;
2927  }
2928 
2929  /* compute infimum = MIN(entropy(op.inf), entropy(op.sup)) and supremum = MAX(MIN(entropy(op.inf), entropy(op.sup))) */
2930 
2931  /* first, compute the logarithms (roundmode nearest, then nextafter) */
2933  if( operand.inf > 0.0 )
2934  {
2935  loginf = log(operand.inf);
2936  infcand1 = SCIPnextafter(loginf, SCIP_REAL_MAX);
2937  supcand1 = SCIPnextafter(loginf, SCIP_REAL_MIN);
2938  }
2939 
2940  if( operand.sup < infinity )
2941  {
2942  logsup = log(operand.sup);
2943  infcand2 = SCIPnextafter(logsup, SCIP_REAL_MAX);
2944  supcand2 = SCIPnextafter(logsup, SCIP_REAL_MIN);
2945  }
2946 
2947  /* second, multiply with operand.inf/sup using upward rounding
2948  * thus, for infinum, negate after muliplication; for supremum, negate before multiplication
2949  */
2951  if( operand.inf > 0.0 )
2952  {
2953  infcand1 = SCIPnegateReal(operand.inf * infcand1);
2954  supcand1 = SCIPnegateReal(operand.inf) * supcand1;
2955  }
2956  else
2957  {
2958  infcand1 = 0.0;
2959  supcand1 = 0.0;
2960  }
2961 
2962  if( operand.sup < infinity )
2963  {
2964  infcand2 = SCIPnegateReal(operand.sup * infcand2);
2965  supcand2 = SCIPnegateReal(operand.sup) * supcand2;
2966  }
2967  else
2968  {
2969  infcand2 = -infinity;
2970  supcand2 = -infinity;
2971  }
2972 
2973  /* restore original rounding mode (asserted to be "to-nearest" above) */
2975 
2976  inf = MIN(infcand1, infcand2);
2977 
2978  extr = exp(-1.0);
2979  if( operand.inf <= extr && extr <= operand.sup )
2980  {
2981  extr = SCIPnextafter(extr, SCIP_REAL_MAX);
2982  sup = MAX3(supcand1, supcand2, extr);
2983  }
2984  else
2985  sup = MAX(supcand1, supcand2);
2986 
2987  assert(inf <= sup);
2988  SCIPintervalSetBounds(resultant, inf, sup);
2989 }
2990 
2991 /** computes exact upper bound on \f$ a x^2 + b x \f$ for x in [xlb, xub], b an interval, and a scalar
2992  *
2993  * Uses Algorithm 2.2 from Domes and Neumaier: Constraint propagation on quadratic constraints (2008).
2994  */
2996  SCIP_Real infinity, /**< value for infinity */
2997  SCIP_Real a, /**< coefficient of x^2 */
2998  SCIP_INTERVAL b_, /**< coefficient of x */
2999  SCIP_INTERVAL x /**< range of x */
3000  )
3001 {
3002  SCIP_Real b;
3003  SCIP_Real u;
3004 
3005  assert(!SCIPintervalIsEmpty(infinity, x));
3006  assert(b_.inf < infinity);
3007  assert(b_.sup > -infinity);
3008  assert( x.inf < infinity);
3009  assert( x.sup > -infinity);
3010 
3011  /* handle b*x separately */
3012  if( a == 0.0 )
3013  {
3014  if( (b_.inf <= -infinity && x.inf < 0.0 ) ||
3015  ( b_.inf < 0.0 && x.inf <= -infinity) ||
3016  ( b_.sup > 0.0 && x.sup >= infinity) ||
3017  ( b_.sup >= infinity && x.sup > 0.0 ) )
3018  {
3019  u = infinity;
3020  }
3021  else
3022  {
3023  SCIP_ROUNDMODE roundmode;
3024  SCIP_Real cand1, cand2, cand3, cand4;
3025 
3026  roundmode = intervalGetRoundingMode();
3028  cand1 = b_.inf * x.inf;
3029  cand2 = b_.inf * x.sup;
3030  cand3 = b_.sup * x.inf;
3031  cand4 = b_.sup * x.sup;
3032  u = MAX(MAX(cand1, cand2), MAX(cand3, cand4));
3033 
3034  intervalSetRoundingMode(roundmode);
3035  }
3036 
3037  return u;
3038  }
3039 
3040  if( x.sup <= 0.0 )
3041  { /* change sign of x: enclose a*x^2 + [-bub, -blb]*(-x) for (-x) in [-xub, -xlb] */
3042  u = x.sup;
3043  x.sup = -x.inf;
3044  x.inf = -u;
3045  b = -b_.inf;
3046  }
3047  else
3048  {
3049  b = b_.sup;
3050  }
3051 
3052  if( x.inf >= 0.0 )
3053  { /* upper bound for a*x^2 + b*x */
3054  SCIP_ROUNDMODE roundmode;
3055  SCIP_Real s,t;
3056 
3057  if( b >= infinity )
3058  return infinity;
3059 
3060  roundmode = intervalGetRoundingMode();
3062 
3063  u = MAX(x.inf * (a*x.inf + b), x.sup * (a*x.sup + b));
3064  s = b/2;
3065  t = s/negate(a);
3066  if( t > x.inf && negate(2*a)*x.sup > b && s*t > u )
3067  u = s*t;
3068 
3069  intervalSetRoundingMode(roundmode);
3070  return u;
3071  }
3072  else
3073  {
3074  SCIP_INTERVAL xlow = x;
3075  SCIP_Real cand1;
3076  SCIP_Real cand2;
3077  assert(x.inf < 0.0 && x.sup > 0);
3078 
3079  xlow.sup = 0; /* so xlow is lower part of interval */
3080  x.inf = 0; /* so x is upper part of interval now */
3081  cand1 = SCIPintervalQuadUpperBound(infinity, a, b_, xlow);
3082  cand2 = SCIPintervalQuadUpperBound(infinity, a, b_, x);
3083  return MAX(cand1, cand2);
3084  }
3085 }
3086 
3087 /** stores range of quadratic term in resultant
3088  *
3089  * given scalar a and intervals b and x, computes interval for \f$ a x^2 + b x \f$ */
3091  SCIP_Real infinity, /**< value for infinity */
3092  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
3093  SCIP_Real sqrcoeff, /**< coefficient of x^2 */
3094  SCIP_INTERVAL lincoeff, /**< coefficient of x */
3095  SCIP_INTERVAL xrng /**< range of x */
3096  )
3097 {
3098  SCIP_Real tmp;
3099 
3100  if( SCIPintervalIsEmpty(infinity, xrng) )
3101  {
3102  SCIPintervalSetEmpty(resultant);
3103  return;
3104  }
3105  if( sqrcoeff == 0.0 )
3106  {
3107  SCIPintervalMul(infinity, resultant, lincoeff, xrng);
3108  return;
3109  }
3110 
3111  resultant->sup = SCIPintervalQuadUpperBound(infinity, sqrcoeff, lincoeff, xrng);
3112 
3113  tmp = lincoeff.inf;
3114  lincoeff.inf = -lincoeff.sup;
3115  lincoeff.sup = -tmp;
3116  resultant->inf = -SCIPintervalQuadUpperBound(infinity, -sqrcoeff, lincoeff, xrng);
3117 
3118  assert(resultant->sup >= resultant->inf);
3119 }
3120 
3121 /** computes interval with positive solutions of a quadratic equation with interval coefficients
3122  *
3123  * Given intervals a, b, and c, this function computes an interval that contains all positive solutions of \f$ a x^2 + b x \in c\f$ within xbnds.
3124  */
3126  SCIP_Real infinity, /**< value for infinity */
3127  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
3128  SCIP_INTERVAL sqrcoeff, /**< coefficient of x^2 */
3129  SCIP_INTERVAL lincoeff, /**< coefficient of x */
3130  SCIP_INTERVAL rhs, /**< right hand side of equation */
3131  SCIP_INTERVAL xbnds /**< bounds on x */
3132  )
3133 {
3134  assert(resultant != NULL);
3135 
3136  /* find x>=0 s.t. a.inf x^2 + b.inf x <= c.sup -> -a.inf x^2 - b.inf x >= -c.sup */
3137  if( lincoeff.inf <= -infinity || rhs.sup >= infinity || sqrcoeff.inf <= -infinity )
3138  {
3139  resultant->inf = 0.0;
3140  resultant->sup = infinity;
3141  }
3142  else
3143  {
3144  SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(infinity, resultant, -sqrcoeff.inf, -lincoeff.inf, -rhs.sup, xbnds);
3145  SCIPdebugMessage("solve %g*x^2 + %g*x >= %g gives [%.20f, %.20f]\n", -sqrcoeff.inf, -lincoeff.inf, -rhs.sup, resultant->inf, resultant->sup);
3146  }
3147 
3148  /* find x>=0 s.t. a.sup x^2 + b.sup x >= c.inf */
3149  if( lincoeff.sup < infinity && rhs.inf > -infinity && sqrcoeff.sup < infinity )
3150  {
3151  SCIP_INTERVAL res2;
3152  /* coverity[uninit_use_in_call] */
3153  SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(infinity, &res2, sqrcoeff.sup, lincoeff.sup, rhs.inf, xbnds);
3154  SCIPdebugMessage("solve %g*x^2 + %g*x >= %g gives [%.20f, %.20f]\n", sqrcoeff.sup, lincoeff.sup, rhs.inf, res2.inf, res2.sup);
3155  SCIPdebugMessage("intersection of [%.20f, %.20f] and [%.20f, %.20f]", resultant->inf, resultant->sup, res2.inf, res2.sup);
3156  /* intersect both results */
3157  SCIPintervalIntersect(resultant, *resultant, res2);
3158  SCIPdebugPrintf(" gives [%.20f, %.20f]\n", resultant->inf, resultant->sup);
3159  }
3160  /* else res2 = [0, infty] */
3161 
3162  if( resultant->inf >= infinity || resultant->sup <= -infinity )
3163  {
3164  SCIPintervalSetEmpty(resultant);
3165  }
3166 }
3167 
3168 /** computes interval with negative solutions of a quadratic equation with interval coefficients
3169  *
3170  * Given intervals a, b, and c, this function computes an interval that contains all negative solutions of \f$ a x^2 + b x \in c\f$ within xbnds.
3171  */
3173  SCIP_Real infinity, /**< value for infinity */
3174  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
3175  SCIP_INTERVAL sqrcoeff, /**< coefficient of x^2 */
3176  SCIP_INTERVAL lincoeff, /**< coefficient of x */
3177  SCIP_INTERVAL rhs, /**< right hand side of equation */
3178  SCIP_INTERVAL xbnds /**< bounds on x */
3179  )
3180 {
3181  SCIP_Real tmp;
3182 
3183  /* change in variables y = -x, thus get all positive solutions of
3184  * a * y^2 + (-b) * y in c with -xbnds as bounds on y
3185  */
3186 
3187  tmp = lincoeff.inf;
3188  lincoeff.inf = -lincoeff.sup;
3189  lincoeff.sup = -tmp;
3190 
3191  tmp = xbnds.inf;
3192  xbnds.inf = -xbnds.sup;
3193  xbnds.sup = -tmp;
3194 
3195  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, resultant, sqrcoeff, lincoeff, rhs, xbnds);
3196 
3197  tmp = resultant->inf;
3198  resultant->inf = -resultant->sup;
3199  resultant->sup = -tmp;
3200 }
3201 
3202 
3203 /** computes positive solutions of a quadratic equation with scalar coefficients
3204  *
3205  * Givens scalar a, b, and c, this function computes an interval that contains all positive solutions of \f$ a x^2 + b x \geq c\f$ within xbnds.
3206  * Implements Algorithm 3.2 from Domes and Neumaier: Constraint propagation on quadratic constraints (2008).
3207  */
3209  SCIP_Real infinity, /**< value for infinity */
3210  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
3211  SCIP_Real sqrcoeff, /**< coefficient of x^2 */
3212  SCIP_Real lincoeff, /**< coefficient of x */
3213  SCIP_Real rhs, /**< right hand side of equation */
3214  SCIP_INTERVAL xbnds /**< bounds on x */
3215  )
3216 {
3217  SCIP_ROUNDMODE roundmode;
3218  SCIP_Real b;
3219  SCIP_Real delta;
3220  SCIP_Real z;
3221 
3222  assert(resultant != NULL);
3223  assert(sqrcoeff < infinity);
3224  assert(sqrcoeff > -infinity);
3225 
3226  if( sqrcoeff == 0.0 )
3227  {
3228  /* special handling for linear b * x >= c
3229  *
3230  * The non-negative solutions here are:
3231  * b < 0, c <= 0 : [0, c/b]
3232  * b <= 0, c > 0 : empty
3233  * b > 0, c > 0 : [c/b, infty]
3234  * b >= 0, c <= 0 : [0, infty]
3235  *
3236  * The same should have been computed below, but without the sqrcoeff, terms simplify (thus, also less rounding).
3237  */
3238 
3239  if( lincoeff <= 0.0 && rhs > 0.0 )
3240  {
3241  SCIPintervalSetEmpty(resultant);
3242  return;
3243  }
3244 
3245  if( lincoeff >= 0.0 && rhs <= 0.0 )
3246  {
3247  /* [0,infty] cap xbnds */
3248  resultant->inf = MAX(0.0, xbnds.inf);
3249  resultant->sup = xbnds.sup;
3250  return;
3251  }
3252 
3253  roundmode = intervalGetRoundingMode();
3254 
3255  if( lincoeff < 0.0 && rhs <= 0.0 )
3256  {
3257  /* [0,c/b] cap xbnds */
3258  resultant->inf = MAX(0.0, xbnds.inf);
3259 
3261  resultant->sup = rhs / lincoeff;
3262  if( xbnds.sup < resultant->sup )
3263  resultant->sup = xbnds.sup;
3264  }
3265  else
3266  {
3267  assert(lincoeff > 0.0);
3268  assert(rhs > 0.0);
3269 
3270  /* [c/b, infty] cap xbnds */
3271 
3273  resultant->inf = rhs / lincoeff;
3274  if( resultant->inf < xbnds.inf )
3275  resultant->inf = xbnds.inf;
3276 
3277  resultant->sup = xbnds.sup;
3278  }
3279 
3280  intervalSetRoundingMode(roundmode);
3281 
3282  return;
3283  }
3284 
3285  resultant->inf = 0.0;
3286  resultant->sup = infinity;
3287 
3288  roundmode = intervalGetRoundingMode();
3289 
3290  /* this should actually be round_upwards, but unless lincoeff is min_double,
3291  * there shouldn't be any rounding happening when dividing by 2, i.e., shifting exponent,
3292  * so it is ok to not change the rounding mode here
3293  */
3294  b = lincoeff / 2.0;
3295 
3296  if( lincoeff >= 0.0 )
3297  { /* b >= 0.0 */
3298  if( rhs > 0.0 )
3299  { /* b >= 0.0 and c > 0.0 */
3301  delta = b*b + sqrcoeff*rhs;
3302  if( delta < 0.0 )
3303  {
3304  SCIPintervalSetEmpty(resultant);
3305  }
3306  else
3307  {
3309  z = SCIPnextafter(sqrt(delta), SCIP_REAL_MAX);
3311  z += b;
3312  resultant->inf = negate(negate(rhs)/z);
3313  if( sqrcoeff < 0.0 )
3314  resultant->sup = z / negate(sqrcoeff);
3315  }
3316  }
3317  else
3318  { /* b >= 0.0 and c <= 0.0 */
3319  if( sqrcoeff < 0.0 )
3320  {
3322  delta = b*b + sqrcoeff*rhs;
3324  z = SCIPnextafter(sqrt(delta), SCIP_REAL_MAX);
3326  z += b;
3327  resultant->sup = z / negate(sqrcoeff);
3328  }
3329  }
3330  }
3331  else
3332  { /* b < 0.0 */
3333  if( rhs > 0.0 )
3334  { /* b < 0.0 and c > 0.0 */
3335  if( sqrcoeff > 0.0 )
3336  {
3338  delta = b*b + sqrcoeff*rhs;
3340  z = SCIPnextafter(sqrt(delta), SCIP_REAL_MIN);
3342  z += negate(b);
3343  resultant->inf = z / sqrcoeff;
3344  }
3345  else
3346  {
3347  SCIPintervalSetEmpty(resultant);
3348  }
3349  }
3350  else
3351  { /* b < 0.0 and c <= 0.0 */
3353  delta = b*b + sqrcoeff * rhs;
3354  if( delta >= 0.0 )
3355  {
3356  /* let resultant = [0,-c/z] for now */
3358  z = SCIPnextafter(sqrt(delta), SCIP_REAL_MIN);
3359  /* continue with downward rounding, because we want z (>= 0) to be small,
3360  * because -rhs/z needs to be large (-rhs >= 0)
3361  */
3363  z += negate(b);
3364  /* also now compute rhs/z with downward rounding, so that -(rhs/z) becomes large */
3365  resultant->sup = negate(rhs/z);
3366 
3367  if( sqrcoeff > 0.0 )
3368  {
3369  /* for a > 0, the result is [0,-c/z] \vee [z/a,infinity]
3370  * currently, resultant = [0,-c/z]
3371  */
3372  SCIP_Real zdiva;
3373 
3374  zdiva = z/sqrcoeff;
3375 
3376  if( xbnds.sup < zdiva )
3377  {
3378  /* after intersecting with xbnds, result is [0,-c/z], so we are done */
3379  }
3380  else if( xbnds.inf > resultant->sup )
3381  {
3382  /* after intersecting with xbnds, result is [z/a,infinity] */
3383  resultant->inf = zdiva;
3384  resultant->sup = infinity;
3385  }
3386  else
3387  {
3388  /* after intersecting with xbnds we can neither exclude [0,-c/z] nor [z/a,infinity],
3389  * so put resultant = [0,infinity] (intersection with xbnds happens below)
3390  * @todo we could create a hole here
3391  */
3392  resultant->sup = infinity;
3393  }
3394  }
3395  else
3396  {
3397  /* for a < 0, the result is [0,-c/z], so we are done */
3398  }
3399  }
3400  }
3401  }
3402 
3403  SCIPintervalIntersect(resultant, *resultant, xbnds);
3404 
3405  intervalSetRoundingMode(roundmode);
3406 }
3407 
3408 /** solves a quadratic equation with interval coefficients
3409  *
3410  * Given intervals a, b and c, this function computes an interval that contains all solutions of \f$ a x^2 + b x \in c\f$ within xbnds.
3411  */
3413  SCIP_Real infinity, /**< value for infinity */
3414  SCIP_INTERVAL* resultant, /**< resultant interval of operation */
3415  SCIP_INTERVAL sqrcoeff, /**< coefficient of x^2 */
3416  SCIP_INTERVAL lincoeff, /**< coefficient of x */
3417  SCIP_INTERVAL rhs, /**< right hand side of equation */
3418  SCIP_INTERVAL xbnds /**< bounds on x */
3419  )
3420 {
3421  SCIP_INTERVAL xpos;
3422  SCIP_INTERVAL xneg;
3423 
3424  assert(resultant != NULL);
3425  assert(!SCIPintervalIsEmpty(infinity, sqrcoeff));
3426  assert(!SCIPintervalIsEmpty(infinity, lincoeff));
3427  assert(!SCIPintervalIsEmpty(infinity, rhs));
3428 
3429  /* special handling for lincoeff * x = rhs without 0 in lincoeff
3430  * rhs/lincoeff gives a good interval that we just have to intersect with xbnds
3431  * the code below would also work, but uses many more case distinctions to get to a result that should be the same (though epsilon differences can sometimes be observed)
3432  */
3433  if( sqrcoeff.inf == 0.0 && sqrcoeff.sup == 0.0 && (lincoeff.inf > 0.0 || lincoeff.sup < 0.0) )
3434  {
3435  SCIPintervalDiv(infinity, resultant, rhs, lincoeff);
3436  SCIPintervalIntersect(resultant, *resultant, xbnds);
3437  SCIPdebugMessage("solving [%g,%g]*x = [%g,%g] for x in [%g,%g] gives [%g,%g]\n", lincoeff.inf, lincoeff.sup, rhs.inf, rhs.sup, xbnds.inf, xbnds.sup, resultant->inf, resultant->sup);
3438  return;
3439  }
3440 
3441  SCIPdebugMessage("solving [%g,%g]*x^2 + [%g,%g]*x = [%g,%g] for x in [%g,%g]\n", sqrcoeff.inf, sqrcoeff.sup, lincoeff.inf, lincoeff.sup, rhs.inf, rhs.sup, xbnds.inf, xbnds.sup);
3442 
3443  /* find all x>=0 such that a*x^2+b*x = c */
3444  if( xbnds.sup >= 0 )
3445  {
3446  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, &xpos, sqrcoeff, lincoeff, rhs, xbnds);
3447  SCIPdebugMessage(" solutions of [%g,%g]*x^2 + [%g,%g]*x in [%g,%g] for x in [%g,%g] are [%.15g,%.15g]\n",
3448  sqrcoeff.inf, sqrcoeff.sup, lincoeff.inf, lincoeff.sup, rhs.inf, rhs.sup, MAX(xbnds.inf, 0.0), xbnds.sup, xpos.inf, xpos.sup);
3449  }
3450  else
3451  {
3452  SCIPintervalSetEmpty(&xpos);
3453  }
3454 
3455  /* find all x<=0 such that a*x^2-b*x = c */
3456  if( xbnds.inf <= 0.0 )
3457  {
3458  SCIPintervalSolveUnivariateQuadExpressionNegative(infinity, &xneg, sqrcoeff, lincoeff, rhs, xbnds);
3459  SCIPdebugMessage(" solutions of [%g,%g]*x^2 + [%g,%g]*x in [%g,%g] for x in [%g,%g] are [%g,%g]\n",
3460  sqrcoeff.inf, sqrcoeff.sup, lincoeff.inf, lincoeff.sup, rhs.inf, rhs.sup, xbnds.inf, MIN(xbnds.sup, 0.0), xneg.inf, xneg.sup);
3461  }
3462  else
3463  {
3464  SCIPintervalSetEmpty(&xneg);
3465  }
3466 
3467  SCIPintervalUnify(resultant, xpos, xneg);
3468  SCIPdebugMessage(" unify gives [%g,%g]\n", SCIPintervalGetInf(*resultant), SCIPintervalGetSup(*resultant));
3469 }
3470 
3471 /** stores range of bivariate quadratic term in resultant
3472  *
3473  * Given scalars \f$a_x\f$, \f$a_y\f$, \f$a_{xy}\f$, \f$b_x\f$, and \f$b_y\f$ and intervals for \f$x\f$ and \f$y\f$,
3474  * computes interval for \f$ a_x x^2 + a_y y^2 + a_{xy} x y + b_x x + b_y y \f$.
3475  *
3476  * \attention The operations are not applied rounding-safe here!
3477  */
3479  SCIP_Real infinity, /**< value for infinity in interval arithmetics */
3480  SCIP_INTERVAL* resultant, /**< buffer where to store result of operation */
3481  SCIP_Real ax, /**< square coefficient of x */
3482  SCIP_Real ay, /**< square coefficient of y */
3483  SCIP_Real axy, /**< bilinear coefficients */
3484  SCIP_Real bx, /**< linear coefficient of x */
3485  SCIP_Real by, /**< linear coefficient of y */
3486  SCIP_INTERVAL xbnds, /**< bounds on x */
3487  SCIP_INTERVAL ybnds /**< bounds on y */
3488  )
3489 {
3490  /* we use double double precision and finally widen the computed range by 1e-8% to compensate for not computing rounding-safe here */
3491  SCIP_Real minval;
3492  SCIP_Real maxval;
3493  SCIP_Real val;
3494  SCIP_Real x;
3495  SCIP_Real y;
3496  SCIP_Real denom;
3497 
3498  assert(resultant != NULL);
3499  assert(xbnds.inf <= xbnds.sup);
3500  assert(ybnds.inf <= ybnds.sup);
3501 
3502  /* if we are separable, then fall back to use SCIPintervalQuad two times and add */
3503  if( axy == 0.0 )
3504  {
3505  SCIP_INTERVAL tmp;
3506 
3507  SCIPintervalSet(&tmp, bx);
3508  SCIPintervalQuad(infinity, resultant, ax, tmp, xbnds);
3509 
3510  SCIPintervalSet(&tmp, by);
3511  SCIPintervalQuad(infinity, &tmp, ay, tmp, ybnds);
3512 
3513  SCIPintervalAdd(infinity, resultant, *resultant, tmp);
3514 
3515  return;
3516  }
3517 
3518  SCIPintervalSet(resultant, 0.0);
3519 
3520  minval = infinity;
3521  maxval = -infinity;
3522 
3523  /* check minima/maxima of expression */
3524  denom = 4.0 * ax * ay - axy * axy;
3525  if( REALABS(denom) > 1e-9 )
3526  {
3527  x = (axy * by - 2.0 * ay * bx) / denom;
3528  y = (axy * bx - 2.0 * ax * by) / denom;
3529  if( xbnds.inf <= x && x <= xbnds.sup && ybnds.inf <= y && y <= ybnds.sup )
3530  {
3531  val = (axy * bx * by - ay * bx * bx - ax * by * by) / denom;
3532  minval = MIN(val, minval);
3533  maxval = MAX(val, maxval);
3534  }
3535  }
3536  else if( REALABS(2.0 * ay * bx - axy * by) <= 1e-9 )
3537  {
3538  /* The whole line (x, -bx/axy - (axy/2ay) x) defines an extreme point with value -ay bx^2 / axy^2
3539  * If x is unbounded, then there is an (x,y) with y in ybnds where the extreme value is assumed.
3540  * 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.
3541  */
3542  if( xbnds.inf <= -infinity && xbnds.sup >= infinity )
3543  {
3544  val = -ay * bx * bx / (axy * axy);
3545  minval = MIN(val, minval);
3546  maxval = MAX(val, maxval);
3547  }
3548  }
3549 
3550  /* check boundary of box xbnds x ybnds */
3551 
3552  if( xbnds.inf <= -infinity )
3553  {
3554  /* check value for x -> -infinity */
3555  if( ax > 0.0 )
3556  maxval = infinity;
3557  else if( ax < 0.0 )
3558  minval = -infinity;
3559  else if( ax == 0.0 )
3560  {
3561  /* bivar(x,y) tends to -(bx+axy y) * infinity */
3562 
3563  if( ybnds.inf <= -infinity )
3564  val = (axy < 0.0 ? -infinity : infinity);
3565  else if( bx + axy * ybnds.inf < 0.0 )
3566  val = infinity;
3567  else
3568  val = -infinity;
3569  minval = MIN(val, minval);
3570  maxval = MAX(val, maxval);
3571 
3572  if( ybnds.sup >= infinity )
3573  val = (axy < 0.0 ? infinity : -infinity);
3574  else if( bx + axy * ybnds.sup < 0.0 )
3575  val = infinity;
3576  else
3577  val = -infinity;
3578  minval = MIN(val, minval);
3579  maxval = MAX(val, maxval);
3580  }
3581  }
3582  else
3583  {
3584  /* get range of bivar(xbnds.inf, y) for y in ybnds */
3585  SCIP_INTERVAL tmp;
3586  SCIP_INTERVAL ycoef;
3587 
3588  SCIPintervalSet(&ycoef, axy * xbnds.inf + by);
3589  SCIPintervalQuad(infinity, &tmp, ay, ycoef, ybnds);
3590  SCIPintervalAddScalar(infinity, &tmp, tmp, (SCIP_Real)(ax * xbnds.inf * xbnds.inf + bx * xbnds.inf));
3591  minval = MIN(tmp.inf, minval);
3592  maxval = MAX(tmp.sup, maxval);
3593  }
3594 
3595  if( xbnds.sup >= infinity )
3596  {
3597  /* check value for x -> infinity */
3598  if( ax > 0.0 )
3599  maxval = infinity;
3600  else if( ax < 0.0 )
3601  minval = -infinity;
3602  else if( ax == 0.0 )
3603  {
3604  /* bivar(x,y) tends to (bx+axy y) * infinity */
3605 
3606  if( ybnds.inf <= -infinity )
3607  val = (axy > 0.0 ? -infinity : infinity);
3608  else if( bx + axy * ybnds.inf > 0.0 )
3609  val = infinity;
3610  else
3611  val = -infinity;
3612  minval = MIN(val, minval);
3613  maxval = MAX(val, maxval);
3614 
3615  if( ybnds.sup >= infinity )
3616  val = (axy > 0.0 ? infinity : -infinity);
3617  else if( bx + axy * ybnds.sup > 0.0 )
3618  val = infinity;
3619  else
3620  val = -infinity;
3621  minval = MIN(val, minval);
3622  maxval = MAX(val, maxval);
3623  }
3624  }
3625  else
3626  {
3627  /* get range of bivar(xbnds.sup, y) for y in ybnds */
3628  SCIP_INTERVAL tmp;
3629  SCIP_INTERVAL ycoef;
3630 
3631  SCIPintervalSet(&ycoef, axy * xbnds.sup + by);
3632  SCIPintervalQuad(infinity, &tmp, ay, ycoef, ybnds);
3633  SCIPintervalAddScalar(infinity, &tmp, tmp, (SCIP_Real)(ax * xbnds.sup * xbnds.sup + bx * xbnds.sup));
3634  minval = MIN(tmp.inf, minval);
3635  maxval = MAX(tmp.sup, maxval);
3636  }
3637 
3638  if( ybnds.inf <= -infinity )
3639  {
3640  /* check value for y -> -infinity */
3641  if( ay > 0.0 )
3642  maxval = infinity;
3643  else if( ay < 0.0 )
3644  minval = -infinity;
3645  else if( ay == 0.0 )
3646  {
3647  /* bivar(x,y) tends to -(by+axy x) * infinity */
3648 
3649  if( xbnds.inf <= -infinity )
3650  val = (axy < 0.0 ? -infinity : infinity);
3651  else if( by + axy * xbnds.inf < 0.0 )
3652  val = infinity;
3653  else
3654  val = -infinity;
3655  minval = MIN(val, minval);
3656  maxval = MAX(val, maxval);
3657 
3658  if( xbnds.sup >= infinity )
3659  val = (axy < 0.0 ? infinity : -infinity);
3660  else if( by + axy * xbnds.sup < 0.0 )
3661  val = infinity;
3662  else
3663  val = -infinity;
3664  minval = MIN(val, minval);
3665  maxval = MAX(val, maxval);
3666  }
3667  }
3668  else
3669  {
3670  /* get range of bivar(x, ybnds.inf) for x in xbnds */
3671  SCIP_INTERVAL tmp;
3672  SCIP_INTERVAL xcoef;
3673 
3674  SCIPintervalSet(&xcoef, axy * ybnds.inf + bx);
3675  SCIPintervalQuad(infinity, &tmp, ax, xcoef, xbnds);
3676  SCIPintervalAddScalar(infinity, &tmp, tmp, ay * ybnds.inf * ybnds.inf + by * ybnds.inf);
3677  minval = MIN(tmp.inf, minval);
3678  maxval = MAX(tmp.sup, maxval);
3679  }
3680 
3681  if( ybnds.sup >= infinity )
3682  {
3683  /* check value for y -> infinity */
3684  if( ay > 0.0 )
3685  maxval = infinity;
3686  else if( ay < 0.0 )
3687  minval = -infinity;
3688  else if( ay == 0.0 )
3689  {
3690  /* bivar(x,y) tends to (by+axy x) * infinity */
3691 
3692  if( xbnds.inf <= -infinity )
3693  val = (axy > 0.0 ? -infinity : infinity);
3694  else if( by + axy * xbnds.inf > 0.0 )
3695  val = infinity;
3696  else
3697  val = -infinity;
3698  minval = MIN(val, minval);
3699  maxval = MAX(val, maxval);
3700 
3701  if( xbnds.sup >= infinity )
3702  val = (axy > 0.0 ? infinity : -infinity);
3703  else if( by + axy * xbnds.sup > 0.0 )
3704  val = infinity;
3705  else
3706  val = -infinity;
3707  minval = MIN(val, minval);
3708  maxval = MAX(val, maxval);
3709  }
3710  }
3711  else
3712  {
3713  /* get range of bivar(x, ybnds.sup) for x in xbnds */
3714  SCIP_INTERVAL tmp;
3715  SCIP_INTERVAL xcoef;
3716 
3717  SCIPintervalSet(&xcoef, axy * ybnds.sup + bx);
3718  SCIPintervalQuad(infinity, &tmp, ax, xcoef, xbnds);
3719  SCIPintervalAddScalar(infinity, &tmp, tmp, (SCIP_Real)(ay * ybnds.sup * ybnds.sup + by * ybnds.sup));
3720  minval = MIN(tmp.inf, minval);
3721  maxval = MAX(tmp.sup, maxval);
3722  }
3723 
3724  minval -= 1e-10 * REALABS(minval);
3725  maxval += 1e-10 * REALABS(maxval);
3726  SCIPintervalSetBounds(resultant, (SCIP_Real)minval, (SCIP_Real)maxval);
3727 
3728  SCIPdebugMessage("range for %gx^2 + %gy^2 + %gxy + %gx + %gy = [%g, %g] for x = [%g, %g], y=[%g, %g]\n",
3729  ax, ay, axy, bx, by, minval, maxval, xbnds.inf, xbnds.sup, ybnds.inf, ybnds.sup);
3730 }
3731 
3732 /** solves a bivariate quadratic equation for the first variable
3733  *
3734  * Given scalars \f$a_x\f$, \f$a_y\f$, \f$a_{xy}\f$, \f$b_x\f$ and \f$b_y\f$, and intervals for \f$x\f$, \f$y\f$, and rhs,
3735  * 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$.
3736  *
3737  * \attention the operations are not applied rounding-safe here
3738  */
3740  SCIP_Real infinity, /**< value for infinity in interval arithmetics */
3741  SCIP_INTERVAL* resultant, /**< buffer where to store result of operation */
3742  SCIP_Real ax, /**< square coefficient of x */
3743  SCIP_Real ay, /**< square coefficient of y */
3744  SCIP_Real axy, /**< bilinear coefficients */
3745  SCIP_Real bx, /**< linear coefficient of x */
3746  SCIP_Real by, /**< linear coefficient of y */
3747  SCIP_INTERVAL rhs, /**< right-hand-side of equation */
3748  SCIP_INTERVAL xbnds, /**< bounds on x */
3749  SCIP_INTERVAL ybnds /**< bounds on y */
3750  )
3751 {
3752  /* we use double double precision and finally widen the computed range by 1e-8% to compensate for not computing rounding-safe here */
3753  SCIP_Real val;
3754 
3755  assert(resultant != NULL);
3756 
3757  if( axy == 0.0 )
3758  {
3759  /* if axy == 0, fall back to SCIPintervalSolveUnivariateQuadExpression */
3760  SCIP_INTERVAL ytermrng;
3761  SCIP_INTERVAL sqrcoef;
3762  SCIP_INTERVAL lincoef;
3763  SCIP_INTERVAL pos;
3764  SCIP_INTERVAL neg;
3765 
3766  SCIPintervalSet(&lincoef, by);
3767  SCIPintervalQuad(infinity, &ytermrng, ay, lincoef, ybnds);
3768  SCIPintervalSub(infinity, &rhs, rhs, ytermrng);
3769 
3770  SCIPintervalSet(&sqrcoef, ax);
3771 
3772  /* get positive solutions, if of interest */
3773  if( xbnds.sup >= 0.0 )
3774  {
3775  SCIPintervalSet(&lincoef, bx);
3776  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, &pos, sqrcoef, lincoef, rhs, xbnds);
3777  }
3778  else
3779  SCIPintervalSetEmpty(&pos);
3780 
3781  /* get negative solutions, if of interest */
3782  if( xbnds.inf < 0.0 )
3783  {
3784  SCIP_INTERVAL xbndsneg;
3785  SCIPintervalSet(&lincoef, -bx);
3786  SCIPintervalSetBounds(&xbndsneg, -xbnds.sup, -xbnds.inf);
3787  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, &neg, sqrcoef, lincoef, rhs, xbndsneg);
3788  if( !SCIPintervalIsEmpty(infinity, neg) )
3789  SCIPintervalSetBounds(&neg, -neg.sup, -neg.inf);
3790  }
3791  else
3792  SCIPintervalSetEmpty(&neg);
3793 
3794  SCIPintervalUnify(resultant, pos, neg);
3795 
3796  return;
3797  }
3798 
3799  if( ybnds.inf <= -infinity || ybnds.sup >= infinity )
3800  {
3801  /* the code below is buggy if y is unbounded, see #2250
3802  * fall back to univariate case by solving a_x x^2 + b_x x + a_y y^2 + (a_xy xbnds + b_y) y in rhs
3803  */
3804  SCIP_INTERVAL ax_;
3805  SCIP_INTERVAL bx_;
3806  SCIP_INTERVAL ycoef;
3807  SCIP_INTERVAL ytermbounds;
3808 
3809  *resultant = xbnds;
3810 
3811  /* nothing we can do here if x is unbounded (we have a_xy != 0 here) */
3812  if( xbnds.inf <= -infinity && xbnds.sup >= infinity )
3813  return;
3814 
3815  /* ycoef = axy xbnds + by */
3816  SCIPintervalMulScalar(infinity, &ycoef, xbnds, axy);
3817  SCIPintervalAddScalar(infinity, &ycoef, ycoef, by);
3818 
3819  /* get bounds on ay y^2 + (axy xbnds + by) y */
3820  SCIPintervalQuad(infinity, &ytermbounds, ay, ycoef, ybnds);
3821 
3822  /* now solve ax x^2 + bx x in rhs - ytermbounds */
3823  SCIPintervalSet(&ax_, ax);
3824  SCIPintervalSet(&bx_, bx);
3825  SCIPintervalSub(infinity, &rhs, rhs, ytermbounds);
3826  SCIPintervalSolveUnivariateQuadExpression(infinity, resultant, ax_, bx_, rhs, xbnds);
3827 
3828  return;
3829  }
3830 
3831  if( ax < 0.0 )
3832  {
3833  SCIP_Real tmp;
3834  tmp = rhs.inf;
3835  rhs.inf = -rhs.sup;
3836  rhs.sup = -tmp;
3837 
3838  SCIPintervalSolveBivariateQuadExpressionAllScalar(infinity, resultant, -ax, -ay, -axy, -bx, -by, rhs, xbnds, ybnds);
3839  return;
3840  }
3841  assert(ax >= 0.0);
3842 
3843  *resultant = xbnds;
3844 
3845  if( ax > 0.0 )
3846  {
3847  SCIP_Real sqrtax;
3848  SCIP_Real minvalleft;
3849  SCIP_Real maxvalleft;
3850  SCIP_Real minvalright;
3851  SCIP_Real maxvalright;
3852  SCIP_Real ymin;
3853  SCIP_Real rcoef_y;
3854  SCIP_Real rcoef_yy;
3855  SCIP_Real rcoef_const;
3856 
3857  sqrtax = sqrt(ax);
3858 
3859  /* rewrite equation as (sqrt(ax)x + b(y))^2 \in r(rhs,y), where
3860  * b(y) = (bx + axy y)/(2sqrt(ax)), r(rhs,y) = rhs - ay y^2 - by y + b(y)^2
3861  *
3862  * -> r(rhs,y) = bx^2/(4ax) + rhs + (axy bx/(2ax) - by)*y + (axy^2/(4ax) - ay)*y^2
3863  */
3864  rcoef_y = axy * bx / (2.0*ax) - by;
3865  rcoef_yy = axy * axy / (4.0*ax) - ay;
3866  rcoef_const = bx * bx / (4.0*ax);
3867 
3868 #define CALCB(y) ((bx + axy * (y)) / (2.0 * sqrtax))
3869 #define CALCR(c,y) (rcoef_const + (c) + (rcoef_y + rcoef_yy * (y)) * (y))
3870 
3871  /* check whether r(rhs,y) is always negative */
3872  if( rhs.sup < infinity )
3873  {
3874  SCIP_INTERVAL ycoef;
3875  SCIP_Real ub;
3876 
3877  SCIPintervalSet(&ycoef, (SCIP_Real)rcoef_y);
3878  ub = (SCIP_Real)(SCIPintervalQuadUpperBound(infinity, (SCIP_Real)rcoef_yy, ycoef, ybnds) + rhs.sup + rcoef_const);
3879 
3880  if( EPSN(ub, 1e-9) )
3881  {
3882  SCIPintervalSetEmpty(resultant);
3883  return;
3884  }
3885  else if( ub < 0.0 )
3886  {
3887  /* it looks like there will be no solution (rhs < 0), but we are very close and above operations did not take care of careful rounding
3888  * thus, we relax rhs a be feasible a bit (-ub would be sufficient, but that would put us exactly onto the boundary)
3889  * see also #1861
3890  */
3891  rhs.sup += -2.0*ub;
3892  }
3893  }
3894 
3895  /* we have sqrt(ax)x \in (-sqrt(r(rhs,y))-b(y)) \cup (sqrt(r(rhs,y))-b(y))
3896  * compute minima and maxima of both functions such that
3897  *
3898  * [minvalleft, maxvalleft ] = -sqrt(r(rhs,y))-b(y)
3899  * [minvalright, maxvalright] = sqrt(r(rhs,y))-b(y)
3900  */
3901 
3902  minvalleft = infinity;
3903  maxvalleft = -infinity;
3904  minvalright = infinity;
3905  maxvalright = -infinity;
3906 
3907  if( rhs.sup >= infinity )
3908  {
3909  /* we can't do much if rhs.sup is infinite
3910  * but we may do a bit of xbnds isn't too huge and rhs.inf > -infinity
3911  */
3912  minvalleft = -infinity;
3913  maxvalright = infinity;
3914  }
3915 
3916  /* evaluate at lower bound of y, as long as r(rhs,ylb) > 0 */
3917  if( ybnds.inf <= -infinity )
3918  {
3919  /* check limit of +/-sqrt(r(rhs,y))-b(y) for y -> -infty */
3920  if( !EPSZ(ay, 1e-9) && axy * axy >= 4.0 * ax * ay )
3921  {
3922  if( axy < 0.0 )
3923  {
3924  minvalleft = -infinity;
3925 
3926  if( ay > 0.0 )
3927  minvalright = -infinity;
3928  else
3929  maxvalright = infinity;
3930  }
3931  else
3932  {
3933  maxvalright = infinity;
3934 
3935  if( ay > 0.0 )
3936  maxvalleft = infinity;
3937  else
3938  minvalleft = -infinity;
3939  }
3940  }
3941  else if( !EPSZ(ay, 1e-9) )
3942  {
3943  /* here axy * axy < 4 * ax * ay, so need to check for zeros of r(rhs,y), which is done below */
3944  }
3945  else
3946  {
3947  /* here ay = 0.0, which gives a limit of -by/2 for -sqrt(r(rhs,y))-b(y) */
3948  minvalleft = -by / 2.0;
3949  maxvalleft = -by / 2.0;
3950  /* here ay = 0.0, which gives a limit of +infinity for sqrt(r(rhs,y))-b(y) */
3951  maxvalright = infinity;
3952  }
3953  }
3954  else
3955  {
3956  SCIP_Real b;
3957  SCIP_Real c;
3958 
3959  b = CALCB(ybnds.inf);
3960 
3961  if( rhs.sup < infinity )
3962  {
3963  c = CALCR(rhs.sup, ybnds.inf);
3964 
3965  if( c > 0.0 )
3966  {
3967  SCIP_Real sqrtc;
3968 
3969  sqrtc = sqrt(c);
3970  minvalleft = MIN(-sqrtc - b, minvalleft);
3971  maxvalright = MAX( sqrtc - b, maxvalright);
3972  }
3973  }
3974 
3975  if( rhs.inf > -infinity )
3976  {
3977  c = CALCR(rhs.inf, ybnds.inf);
3978 
3979  if( c > 0.0 )
3980  {
3981  SCIP_Real sqrtc;
3982 
3983  sqrtc = sqrt(c);
3984  maxvalleft = MAX(-sqrtc - b, maxvalleft);
3985  minvalright = MIN( sqrtc - b, minvalright);
3986  }
3987  }
3988  }
3989 
3990  /* evaluate at upper bound of y, as long as r(rhs, yub) > 0 */
3991  if( ybnds.sup >= infinity )
3992  {
3993  /* check limit of +/-sqrt(r(rhs,y))-b(y) for y -> +infty */
3994  if( !EPSZ(ay, 1e-9) && axy * axy >= 4.0 * ax * ay )
3995  {
3996  if( axy > 0.0 )
3997  {
3998  minvalleft = -infinity;
3999 
4000  if( ay > 0.0 )
4001  minvalright = -infinity;
4002  else
4003  maxvalright = infinity;
4004  }
4005  else
4006  {
4007  maxvalright = infinity;
4008 
4009  if( ay > 0.0 )
4010  maxvalleft = infinity;
4011  else
4012  minvalleft = -infinity;
4013  }
4014  }
4015  else if( !EPSZ(ay, 1e-9) )
4016  {
4017  /* here axy * axy < 4 * ax * ay, so need to check for zeros of r(rhs,y), which will happen below */
4018  }
4019  else
4020  {
4021  /* here ay = 0.0, which gives a limit of -infinity for -sqrt(r(rhs,y))-b(y) */
4022  minvalleft = -infinity;
4023  /* here ay = 0.0, which gives a limit of -by/2 for sqrt(r(rhs,y))-b(y) */
4024  minvalright = MIN(minvalright, -by / 2.0);
4025  maxvalright = MAX(maxvalright, -by / 2.0);
4026  }
4027  }
4028  else
4029  {
4030  SCIP_Real b;
4031  SCIP_Real c;
4032 
4033  b = CALCB(ybnds.sup);
4034 
4035  if( rhs.sup < infinity )
4036  {
4037  c = CALCR(rhs.sup, ybnds.sup);
4038 
4039  if( c > 0.0 )
4040  {
4041  SCIP_Real sqrtc;
4042 
4043  sqrtc = sqrt(c);
4044  minvalleft = MIN(-sqrtc - b, minvalleft);
4045  maxvalright = MAX( sqrtc - b, maxvalright);
4046  }
4047  }
4048 
4049  if( rhs.inf > -infinity )
4050  {
4051  c = CALCR(rhs.inf, ybnds.sup);
4052 
4053  if( c > 0.0 )
4054  {
4055  SCIP_Real sqrtc;
4056 
4057  sqrtc = sqrt(c);
4058  maxvalleft = MAX(-sqrtc - b, maxvalleft);
4059  minvalright = MIN( sqrtc - b, minvalright);
4060  }
4061  }
4062  }
4063 
4064  /* evaluate at ymin = y_{_,+}, if inside ybnds
4065  * if ay = 0 or 2ay*bx == axy*by, then there is no ymin */
4066  if( !EPSZ(ay, 1e-9) )
4067  {
4068  if( REALABS(axy*axy - 4.0*ax*ay) > 1e-9 )
4069  {
4070  SCIP_Real sqrtterm;
4071 
4072  if( rhs.sup < infinity )
4073  {
4074  sqrtterm = axy * axy * ay * (ay * bx * bx - axy * bx * by + ax * by * by - axy * axy * rhs.sup + 4.0 * ax * ay * rhs.sup);
4075  if( !EPSN(sqrtterm, 1e-9) )
4076  {
4077  sqrtterm = sqrt(MAX(sqrtterm, 0.0));
4078  /* check first candidate for extreme points of +/-sqrt(rhs(r,y))-b(y) */
4079  ymin = axy * ay * bx - 2.0 * ax * ay * by - sqrtterm;
4080  ymin /= ay;
4081  ymin /= 4.0 * ax * ay - axy * axy;
4082 
4083  if( ymin > ybnds.inf && ymin < ybnds.sup )
4084  {
4085  SCIP_Real b;
4086  SCIP_Real c;
4087 
4088  b = CALCB(ymin);
4089  c = CALCR(rhs.sup, ymin);
4090 
4091  if( c > 0.0 )
4092  {
4093  SCIP_Real sqrtc;
4094 
4095  sqrtc = sqrt(c);
4096  minvalleft = MIN(-sqrtc - b, minvalleft);
4097  maxvalright = MAX( sqrtc - b, maxvalright);
4098  }
4099  }
4100 
4101  /* check second candidate for extreme points of +/-sqrt(rhs(r,y))-b(y) */
4102  ymin = axy * ay * bx - 2.0 * ax * ay * by + sqrtterm;
4103  ymin /= ay;
4104  ymin /= 4.0 * ax * ay - axy * axy;
4105 
4106  if( ymin > ybnds.inf && ymin < ybnds.sup )
4107  {
4108  SCIP_Real b;
4109  SCIP_Real c;
4110 
4111  b = CALCB(ymin);
4112  c = CALCR(rhs.sup, ymin);
4113 
4114  if( c > 0.0 )
4115  {
4116  SCIP_Real sqrtc;
4117 
4118  sqrtc = sqrt(c);
4119  minvalleft = MIN(-sqrtc - b, minvalleft);
4120  maxvalright = MAX( sqrtc - b, maxvalright);
4121  }
4122  }
4123  }
4124  }
4125 
4126  if( rhs.inf > -infinity )
4127  {
4128  sqrtterm = axy * axy * ay * (ay * bx * bx - axy * bx * by + ax * by * by - axy * axy * rhs.inf + 4.0 * ax * ay * rhs.inf);
4129  if( !EPSN(sqrtterm, 1e-9) )
4130  {
4131  sqrtterm = sqrt(MAX(sqrtterm, 0.0));
4132  /* check first candidate for extreme points of +/-sqrt(r(rhs,y))-b(y) */
4133  ymin = axy * ay * bx - 2.0 * ax * ay * by - sqrtterm;
4134  ymin /= ay;
4135  ymin /= 4.0 * ax * ay - axy * axy;
4136 
4137  if( ymin > ybnds.inf && ymin < ybnds.sup )
4138  {
4139  SCIP_Real b;
4140  SCIP_Real c;
4141 
4142  b = CALCB(ymin);
4143  c = CALCR(rhs.inf, ymin);
4144 
4145  if( c > 0.0 )
4146  {
4147  SCIP_Real sqrtc;
4148 
4149  sqrtc = sqrt(c);
4150  maxvalleft = MAX(-sqrtc - b, maxvalleft);
4151  minvalright = MIN( sqrtc - b, minvalright);
4152  }
4153  }
4154 
4155  /* check second candidate for extreme points of +/-sqrt(c(y))-b(y) */
4156  ymin = axy * ay * bx - 2.0 * ax * ay * by + sqrtterm;
4157  ymin /= ay;
4158  ymin /= 4.0 * ax * ay - axy * axy;
4159 
4160  if( ymin > ybnds.inf && ymin < ybnds.sup )
4161  {
4162  SCIP_Real b;
4163  SCIP_Real c;
4164 
4165  b = CALCB(ymin);
4166  c = CALCR(rhs.inf, ymin);
4167 
4168  if( c > 0.0 )
4169  {
4170  SCIP_Real sqrtc;
4171 
4172  sqrtc = sqrt(c);
4173  maxvalleft = MAX(-sqrtc - b, maxvalleft);
4174  minvalright = MIN( sqrtc - b, minvalright);
4175  }
4176  }
4177  }
4178  }
4179  }
4180  else if( REALABS(2.0 * ay * bx - axy * by) > 1e-9 )
4181  {
4182  if( rhs.sup < infinity )
4183  {
4184  ymin = - (4.0 * ay * bx * by - axy * by * by + 4.0 * axy * ay * rhs.sup);
4185  ymin /= 4.0 * ay;
4186  ymin /= 2.0 * ay * bx - axy * by;
4187 
4188  if( ymin > ybnds.inf && ymin < ybnds.sup )
4189  {
4190  SCIP_Real b;
4191  SCIP_Real c;
4192 
4193  b = CALCB(ymin);
4194  c = CALCR(rhs.sup, ymin);
4195 
4196  if( c > 0.0 )
4197  {
4198  SCIP_Real sqrtc;
4199 
4200  sqrtc = sqrt(c);
4201  minvalleft = MIN(-sqrtc - b, minvalleft);
4202  maxvalright = MAX( sqrtc - b, maxvalright);
4203  }
4204  }
4205  }
4206 
4207  if( rhs.inf > -infinity )
4208  {
4209  ymin = - (4.0 * ay * bx * by - axy * by * by + 4.0 * axy * ay * rhs.inf);
4210  ymin /= 4.0 * ay;
4211  ymin /= 2.0 * ay * bx - axy * by;
4212 
4213  if( ymin > ybnds.inf && ymin < ybnds.sup )
4214  {
4215  SCIP_Real b;
4216  SCIP_Real c;
4217 
4218  b = CALCB(ymin);
4219  c = CALCR(rhs.inf, ymin);
4220 
4221  if( c > 0.0 )
4222  {
4223  SCIP_Real sqrtc;
4224 
4225  sqrtc = sqrt(c);
4226  maxvalleft = MAX(-sqrtc - b, maxvalleft);
4227  minvalright = MIN( sqrtc - b, minvalright);
4228  }
4229  }
4230  }
4231  }
4232  }
4233 
4234  /* evaluate the case r(rhs,y) = 0, which is to min/max -b(y) w.r.t. r(rhs,y) = 0, y in ybnds
4235  * with the above assignments
4236  * rcoef_y = axy * bx / (2.0*ax) - by;
4237  * rcoef_yy = axy * axy / (4.0*ax) - ay;
4238  * rcoef_const = bx * bx / (4.0*ax);
4239  * we have r(rhs,y) = rhs + rcoef_const + rcoef_y * y + rcoef_yy * y^2
4240  *
4241  * thus, r(rhs,y) = 0 <-> rcoef_y * y + rcoef_yy * y^2 in -rhs - rcoef_const
4242  *
4243  */
4244  {
4245  SCIP_INTERVAL rcoef_yy_int;
4246  SCIP_INTERVAL rcoef_y_int;
4247  SCIP_INTERVAL rhs2;
4248  SCIP_Real b;
4249 
4250  /* setup rcoef_yy, rcoef_y and -rhs-rcoef_const as intervals */
4251  SCIPintervalSet(&rcoef_yy_int, (SCIP_Real)rcoef_yy);
4252  SCIPintervalSet(&rcoef_y_int, (SCIP_Real)rcoef_y);
4253  SCIPintervalSetBounds(&rhs2, (SCIP_Real)(-rhs.sup - rcoef_const), (SCIP_Real)(-rhs.inf - rcoef_const));
4254 
4255  /* first find all y >= 0 such that rcoef_y * y + rcoef_yy * y^2 in -rhs2, if ybnds.sup >= 0.0
4256  * and evaluate -b(y) w.r.t. these values
4257  */
4258  if( ybnds.sup >= 0.0 )
4259  {
4260  SCIP_INTERVAL ypos;
4261 
4262  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, &ypos, rcoef_yy_int, rcoef_y_int, rhs2, ybnds);
4263  if( !SCIPintervalIsEmpty(infinity, ypos) )
4264  {
4265  assert(ypos.inf >= 0.0); /* we computed only positive solutions above */
4266  b = CALCB(ypos.inf);
4267  minvalleft = MIN(minvalleft, -b);
4268  maxvalleft = MAX(maxvalleft, -b);
4269  minvalright = MIN(minvalright, -b);
4270  maxvalright = MAX(maxvalright, -b);
4271 
4272  if( ypos.sup < infinity )
4273  {
4274  b = CALCB(ypos.sup);
4275  minvalleft = MIN(minvalleft, -b);
4276  maxvalleft = MAX(maxvalleft, -b);
4277  minvalright = MIN(minvalright, -b);
4278  maxvalright = MAX(maxvalright, -b);
4279  }
4280  else
4281  {
4282  /* -b(y) = - (bx + axy * y) / (2.0 * sqrt(ax)) -> -sign(axy)*infinity for y -> infinity */
4283  if( axy > 0.0 )
4284  {
4285  minvalleft = -infinity;
4286  minvalright = -infinity;
4287  }
4288  else
4289  {
4290  maxvalleft = infinity;
4291  maxvalright = infinity;
4292  }
4293  }
4294  }
4295  }
4296 
4297  /* next find all y <= 0 such that rcoef_y * y + rcoef_yy * y^2 in -rhs2, if ybnds.inf < 0.0
4298  * and evaluate -b(y) w.r.t. these values
4299  * (the case y fixed to 0 has been handled in the ybnds.sup >= 0 case above)
4300  */
4301  if( ybnds.inf < 0.0 )
4302  {
4303  SCIP_INTERVAL yneg;
4304 
4305  SCIPintervalSolveUnivariateQuadExpressionNegative(infinity, &yneg, rcoef_yy_int, rcoef_y_int, rhs2, ybnds);
4306  if( !SCIPintervalIsEmpty(infinity, yneg) )
4307  {
4308  if( yneg.inf > -infinity )
4309  {
4310  b = CALCB(yneg.inf);
4311  minvalleft = MIN(minvalleft, -b);
4312  maxvalleft = MAX(maxvalleft, -b);
4313  minvalright = MIN(minvalright, -b);
4314  maxvalright = MAX(maxvalright, -b);
4315  }
4316  else
4317  {
4318  /* -b(y) = - (bx + axy * y) / (2.0 * sqrt(ax)) -> sign(axy)*infinity for y -> -infinity */
4319  if( axy > 0.0 )
4320  {
4321  maxvalleft = infinity;
4322  maxvalright = infinity;
4323  }
4324  else
4325  {
4326  minvalleft = -infinity;
4327  minvalright = -infinity;
4328  }
4329  }
4330 
4331  assert(yneg.sup <= 0.0); /* we computed only negative solutions above */
4332  b = CALCB(yneg.sup);
4333  minvalleft = MIN(minvalleft, -b);
4334  maxvalleft = MAX(maxvalleft, -b);
4335  minvalright = MIN(minvalright, -b);
4336  maxvalright = MAX(maxvalright, -b);
4337  }
4338  }
4339  }
4340 
4341  if( rhs.inf > -infinity && xbnds.inf > -infinity && EPSGT(xbnds.inf, maxvalleft / sqrtax, 1e-9) )
4342  {
4343  /* 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)
4344  * 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) */
4345  assert(EPSGE(minvalright, minvalleft, 1e-9)); /* right interval should not be above lower bound of left interval */
4346  if( minvalright > -infinity )
4347  {
4348  assert(minvalright < infinity);
4349  resultant->inf = (SCIP_Real)(minvalright / sqrtax);
4350  }
4351  }
4352  else
4353  {
4354  /* otherwise, tighten lower bound of sqrt(ax)*x to lower bound of -sqrt(r(rhs,y))-b(y) */
4355  if( minvalleft > -infinity )
4356  {
4357  assert(minvalleft < infinity);
4358  resultant->inf = (SCIP_Real)(minvalleft / sqrtax);
4359  }
4360  }
4361 
4362  if( rhs.inf > -infinity && xbnds.sup < infinity && EPSLT(xbnds.sup, minvalright / sqrtax, 1e-9) )
4363  {
4364  /* 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)
4365  * 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) */
4366  assert(EPSLE(maxvalleft, maxvalright, 1e-9)); /* left interval should not be above upper bound of right interval */
4367  if( maxvalleft < infinity )
4368  {
4369  assert(maxvalleft > -infinity);
4370  resultant->sup = (SCIP_Real)(maxvalleft / sqrtax);
4371  }
4372  }
4373  else
4374  {
4375  /* otherwise, tighten upper bound of sqrt(ax)*x to upper bound of sqrt(r(rhs,y))-b(y) */
4376  if( maxvalright < infinity )
4377  {
4378  assert(maxvalright > -infinity);
4379  resultant->sup = (SCIP_Real)(maxvalright / sqrtax);
4380  }
4381  }
4382 
4383  resultant->inf -= 1e-10 * REALABS(resultant->inf);
4384  resultant->sup += 1e-10 * REALABS(resultant->sup);
4385 
4386 #undef CALCB
4387 #undef CALCR
4388  }
4389  else
4390  {
4391  /* case ax == 0 */
4392 
4393  SCIP_Real c;
4394  SCIP_Real d;
4395  SCIP_Real ymin;
4396  SCIP_Real minval;
4397  SCIP_Real maxval;
4398 
4399  /* consider -bx / axy in ybnds, i.e., bx + axy y can be 0 */
4400  if( EPSGE(-bx / axy, ybnds.inf, 1e-9) && EPSLE(-bx / axy, ybnds.sup, 1e-9) )
4401  {
4402  /* write as (bx + axy y) * x \in (c - ay y^2 - by y)
4403  * and estimate bx + axy y and c - ay y^2 - by y by intervals independently
4404  * @todo can we do better, as in the case where bx + axy y is bounded away from 0?
4405  */
4406  SCIP_INTERVAL lincoef;
4407  SCIP_INTERVAL myrhs;
4408  SCIP_INTERVAL tmp;
4409 
4410  if( xbnds.inf < 0.0 && xbnds.sup > 0.0 )
4411  {
4412  /* if (bx + axy y) can be arbitrary small and x be both positive and negative,
4413  * then nothing we can tighten here
4414  */
4415  SCIPintervalSetBounds(resultant, xbnds.inf, xbnds.sup);
4416  return;
4417  }
4418 
4419  /* store interval for (bx + axy y) in lincoef */
4420  SCIPintervalMulScalar(infinity, &lincoef, ybnds, axy);
4421  SCIPintervalAddScalar(infinity, &lincoef, lincoef, bx);
4422 
4423  /* store interval for (c - ay y^2 - by y) in myrhs */
4424  SCIPintervalSet(&tmp, by);
4425  SCIPintervalQuad(infinity, &tmp, ay, tmp, ybnds);
4426  SCIPintervalSub(infinity, &myrhs, rhs, tmp);
4427 
4428  if( lincoef.inf == 0.0 && lincoef.sup == 0.0 )
4429  {
4430  /* equation became 0.0 \in myrhs */
4431  if( myrhs.inf <= 0.0 && myrhs.sup >= 0.0 )
4432  *resultant = xbnds;
4433  else
4434  SCIPintervalSetEmpty(resultant);
4435  }
4436  else if( xbnds.inf >= 0.0 )
4437  {
4438  SCIP_INTERVAL a_;
4439 
4440  /* need only positive solutions */
4441  SCIPintervalSet(&a_, 0.0);
4442  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, resultant, a_, lincoef, myrhs, xbnds);
4443  }
4444  else
4445  {
4446  SCIP_INTERVAL a_;
4447  SCIP_INTERVAL xbndsneg;
4448 
4449  assert(xbnds.sup <= 0.0);
4450 
4451  /* need only negative solutions */
4452  SCIPintervalSet(&a_, 0.0);
4453  SCIPintervalSetBounds(&lincoef, -lincoef.sup, -lincoef.inf);
4454  SCIPintervalSetBounds(&xbndsneg, -xbnds.sup, -xbnds.inf);
4455  SCIPintervalSolveUnivariateQuadExpressionPositive(infinity, resultant, a_, lincoef, myrhs, xbndsneg);
4456  if( !SCIPintervalIsEmpty(infinity, *resultant) )
4457  SCIPintervalSetBounds(resultant, -resultant->sup, -resultant->inf);
4458  }
4459 
4460  return;
4461  }
4462 
4463  minval = infinity;
4464  maxval = -infinity;
4465 
4466  /* compute a lower bound on x */
4467  if( bx + axy * (axy > 0.0 ? ybnds.inf : ybnds.sup) > 0.0 )
4468  c = rhs.inf;
4469  else
4470  c = rhs.sup;
4471 
4472  if( c > -infinity && c < infinity )
4473  {
4474  if( ybnds.inf <= -infinity )
4475  {
4476  /* limit is ay/axy * infinity if ay != 0.0 and -by/axy otherwise */
4477  if( EPSZ(ay, 1e-9) )
4478  minval = -by / axy;
4479  else if( ay * axy < 0.0 )
4480  minval = -infinity;
4481  }
4482  else
4483  {
4484  val = (c - ay * ybnds.inf * ybnds.inf - by * ybnds.inf) / (bx + axy * ybnds.inf);
4485  minval = MIN(val, minval);
4486  }
4487 
4488  if( ybnds.sup >= infinity )
4489  {
4490  /* limit is -ay/axy * infinity if ay != 0.0 and -by/axy otherwise */
4491  if( EPSZ(ay, 1e-9) )
4492  minval = MIN(minval, -by / axy);
4493  else if( ay * axy > 0.0 )
4494  minval = -infinity;
4495  }
4496  else
4497  {
4498  val = (c - ay * ybnds.sup * ybnds.sup - by * ybnds.sup) / (bx + axy * ybnds.sup);
4499  minval = MIN(val, minval);
4500  }
4501 
4502  if( !EPSZ(ay, 1e-9) )
4503  {
4504  d = ay * (ay * bx * bx - axy * (bx * by + axy * c));
4505  if( !EPSN(d, 1e-9) )
4506  {
4507  ymin = -ay * bx + sqrt(MAX(d, 0.0));
4508  ymin /= axy * ay;
4509 
4510  if( ymin > ybnds.inf && ymin < ybnds.sup )
4511  {
4512  assert(bx + axy * ymin != 0.0);
4513 
4514  val = (c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4515  minval = MIN(val, minval);
4516  }
4517 
4518  ymin = -ay * bx - sqrt(MAX(d, 0.0));
4519  ymin /= axy * ay;
4520 
4521  if(ymin > ybnds.inf && ymin < ybnds.sup )
4522  {
4523  assert(bx + axy * ymin != 0.0);
4524 
4525  val = (c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4526  minval = MIN(val, minval);
4527  }
4528  }
4529  }
4530  }
4531  else
4532  {
4533  minval = -infinity;
4534  }
4535 
4536  /* compute an upper bound on x */
4537  if( bx + axy * (axy > 0.0 ? ybnds.inf : ybnds.sup) > 0.0 )
4538  c = rhs.sup;
4539  else
4540  c = rhs.inf;
4541 
4542  if( c > -infinity && c < infinity )
4543  {
4544  if( ybnds.inf <= -infinity )
4545  {
4546  /* limit is ay/axy * infinity if ay != 0.0 and -by/axy otherwise */
4547  if( EPSZ(ay, 1e-9) )
4548  maxval = -by / axy;
4549  else if( ay * axy > 0.0 )
4550  maxval = infinity;
4551  }
4552  else
4553  {
4554  val = (c - ay * ybnds.inf * ybnds.inf - by * ybnds.inf) / (bx + axy * ybnds.inf);
4555  maxval = MAX(val, maxval);
4556  }
4557 
4558  if( ybnds.sup >= infinity )
4559  {
4560  /* limit is -ay/axy * infinity if ay != 0.0 and -by/axy otherwise */
4561  if( EPSZ(ay, 1e-9) )
4562  maxval = MAX(maxval, -by / axy);
4563  else if( ay * axy < 0.0 )
4564  maxval = infinity;
4565  }
4566  else
4567  {
4568  val = (c - ay * ybnds.sup * ybnds.sup - by * ybnds.sup) / (bx + axy * ybnds.sup);
4569  maxval = MAX(val, maxval);
4570  }
4571 
4572  if( !EPSZ(ay, 1e-9) )
4573  {
4574  d = ay * (ay * bx * bx - axy * (bx * by + axy * c));
4575  if( !EPSN(d, 1e-9) )
4576  {
4577  ymin = ay * bx + sqrt(MAX(d, 0.0));
4578  ymin /= axy * ay;
4579 
4580  if( ymin > ybnds.inf && ymin < ybnds.sup )
4581  {
4582  assert(bx + axy * ymin != 0.0); /* the case -bx/axy in ybnds was handled aboved */
4583  val = (c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4584  maxval = MAX(val, maxval);
4585  }
4586 
4587  ymin = ay * bx - sqrt(MAX(d, 0.0));
4588  ymin /= axy * ay;
4589 
4590  if( ymin > ybnds.inf && ymin < ybnds.sup )
4591  {
4592  assert(bx + axy * ymin != 0.0); /* the case -bx/axy in ybnds was handled aboved */
4593  val = (c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4594  maxval = MAX(val, maxval);
4595  }
4596  }
4597  }
4598  }
4599  else
4600  {
4601  maxval = infinity;
4602  }
4603 
4604  if( minval > -infinity )
4605  resultant->inf = minval - 1e-10 * REALABS(minval);
4606  else
4607  resultant->inf = -infinity;
4608  if( maxval < infinity )
4609  resultant->sup = maxval + 1e-10 * REALABS(maxval);
4610  else
4611  resultant->sup = infinity;
4612  SCIPintervalIntersect(resultant, *resultant, xbnds);
4613  }
4614 }
4615 
4616 /** propagates a weighted sum of intervals in a given interval
4617  *
4618  * Given \f$\text{constant} + \sum_i \text{weights}_i \text{operands}_i \in \text{rhs}\f$,
4619  * computes possibly tighter interval for each term.
4620  *
4621  * @attention Valid values are returned in resultants only if any tightening has been found and no empty interval, that is, function returns with non-zero and `*infeasible` = FALSE.
4622  *
4623  * @return Number of terms for which resulting interval is smaller than operand interval.
4624  */
4626  SCIP_Real infinity, /**< value for infinity in interval arithmetics */
4627  int noperands, /**< number of operands (intervals) to propagate */
4628  SCIP_INTERVAL* operands, /**< intervals to propagate */
4629  SCIP_Real* weights, /**< weights of intervals in sum */
4630  SCIP_Real constant, /**< constant in sum */
4631  SCIP_INTERVAL rhs, /**< right-hand-side interval */
4632  SCIP_INTERVAL* resultants, /**< array to store propagated intervals, if any reduction is found at all (check return code and *infeasible) */
4633  SCIP_Bool* infeasible /**< buffer to store if propagation produced empty interval */
4634  )
4635 {
4636  SCIP_ROUNDMODE prevroundmode;
4637  SCIP_INTERVAL childbounds;
4638  SCIP_Real minlinactivity;
4639  SCIP_Real maxlinactivity;
4640  int minlinactivityinf;
4641  int maxlinactivityinf;
4642  int nreductions = 0;
4643  int c;
4644 
4645  assert(noperands > 0);
4646  assert(operands != NULL);
4647  assert(weights != NULL);
4648  assert(resultants != NULL);
4649  assert(infeasible != NULL);
4650 
4651  *infeasible = FALSE;
4652 
4653  /* not possible to conclude finite bounds if the rhs is [-inf,inf] */
4654  if( SCIPintervalIsEntire(infinity, rhs) )
4655  return 0;
4656 
4657  prevroundmode = SCIPintervalGetRoundingMode();
4659 
4660  minlinactivity = constant;
4661  maxlinactivity = -constant; /* use -constant because of the rounding mode */
4662  minlinactivityinf = 0;
4663  maxlinactivityinf = 0;
4664 
4665  SCIPdebugMessage("reverse prop with %d children: %.20g", noperands, constant);
4666 
4667  /* shift coefficients into the intervals of the children (using resultants as working memory)
4668  * compute the min and max activities
4669  */
4670  for( c = 0; c < noperands; ++c )
4671  {
4672  childbounds = operands[c];
4673  SCIPdebugPrintf(" %+.20g*[%.20g,%.20g]", weights[c], childbounds.inf, childbounds.sup);
4674 
4675  if( SCIPintervalIsEmpty(infinity, childbounds) )
4676  {
4677  *infeasible = TRUE;
4678  c = noperands; /* signal for terminate code to not copy operands to resultants because we return *infeasible == TRUE */ /*lint !e850*/
4679  goto TERMINATE;
4680  }
4681 
4682  SCIPintervalMulScalar(infinity, &resultants[c], childbounds, weights[c]);
4683 
4684  if( resultants[c].sup >= infinity )
4685  ++maxlinactivityinf;
4686  else
4687  {
4688  assert(resultants[c].sup > -infinity);
4689  maxlinactivity -= resultants[c].sup;
4690  }
4691 
4692  if( resultants[c].inf <= -infinity )
4693  ++minlinactivityinf;
4694  else
4695  {
4696  assert(resultants[c].inf < infinity);
4697  minlinactivity += resultants[c].inf;
4698  }
4699  }
4700  maxlinactivity = -maxlinactivity; /* correct sign */
4701 
4702  SCIPdebugPrintf(" = [%.20g,%.20g] in rhs = [%.20g,%.20g]\n",
4703  minlinactivityinf ? -infinity : minlinactivity,
4704  maxlinactivityinf ? infinity : maxlinactivity,
4705  rhs.inf, rhs.sup);
4706 
4707  /* if there are too many unbounded bounds, then could only compute infinite bounds for children, so give up */
4708  if( (minlinactivityinf >= 2 || rhs.sup >= infinity) && (maxlinactivityinf >= 2 || rhs.inf <= -infinity) )
4709  {
4710  c = noperands; /* signal for terminate code that it doesn't need to copy operands to resultants because we return nreductions==0 */
4711  goto TERMINATE;
4712  }
4713 
4714  for( c = 0; c < noperands; ++c )
4715  {
4716  /* upper bounds of c_i is
4717  * node->bounds.sup - (minlinactivity - c_i.inf), if c_i.inf > -infinity and minlinactivityinf == 0
4718  * node->bounds.sup - minlinactivity, if c_i.inf == -infinity and minlinactivityinf == 1
4719  */
4720  SCIPintervalSetEntire(infinity, &childbounds);
4721  if( rhs.sup < infinity )
4722  {
4723  /* we are still in downward rounding mode, so negate and negate to get upward rounding */
4724  if( resultants[c].inf <= -infinity && minlinactivityinf <= 1 )
4725  {
4726  assert(minlinactivityinf == 1);
4727  childbounds.sup = SCIPintervalNegateReal(minlinactivity - rhs.sup);
4728  }
4729  else if( minlinactivityinf == 0 )
4730  {
4731  childbounds.sup = SCIPintervalNegateReal(minlinactivity - rhs.sup - resultants[c].inf);
4732  }
4733  }
4734 
4735  /* lower bounds of c_i is
4736  * node->bounds.inf - (maxlinactivity - c_i.sup), if c_i.sup < infinity and maxlinactivityinf == 0
4737  * node->bounds.inf - maxlinactivity, if c_i.sup == infinity and maxlinactivityinf == 1
4738  */
4739  if( rhs.inf > -infinity )
4740  {
4741  if( resultants[c].sup >= infinity && maxlinactivityinf <= 1 )
4742  {
4743  assert(maxlinactivityinf == 1);
4744  childbounds.inf = rhs.inf - maxlinactivity;
4745  }
4746  else if( maxlinactivityinf == 0 )
4747  {
4748  childbounds.inf = rhs.inf - maxlinactivity + resultants[c].sup;
4749  }
4750  }
4751 
4752  SCIPdebugMessage("child %d: %.20g*x in [%.20g,%.20g]", c, weights[c], childbounds.inf, childbounds.sup);
4753 
4754  /* if interval arithmetics works correctly, then childbounds cannot be empty, which is also asserted in SCIPintervalDivScalar below
4755  * however, if running under valgrind, then interval arithmetics does not work correctly and thus one may run into an assert in SCIPintervalDivScalar
4756  * so this check avoids this by declaring the propagation to result in an empty interval (thus only do if asserts are on)
4757  */
4758 #ifndef NDEBUG
4759  if( SCIPintervalIsEmpty(infinity, childbounds) )
4760  {
4761  *infeasible = TRUE;
4762  c = noperands; /*lint !e850*/
4763  goto TERMINATE;
4764  }
4765 #endif
4766 
4767  /* divide by the child coefficient */
4768  SCIPintervalDivScalar(infinity, &childbounds, childbounds, weights[c]);
4769 
4770  SCIPdebugPrintf(" -> x = [%.20g,%.20g]\n", childbounds.inf, childbounds.sup);
4771 
4772  SCIPintervalIntersect(&resultants[c], operands[c], childbounds);
4773  if( SCIPintervalIsEmpty(infinity, resultants[c]) )
4774  {
4775  *infeasible = TRUE;
4776  c = noperands; /*lint !e850*/
4777  goto TERMINATE;
4778  }
4779 
4780  if( resultants[c].inf != operands[c].inf || resultants[c].sup != operands[c].sup )
4781  ++nreductions;
4782  }
4783 
4784 TERMINATE:
4785  SCIPintervalSetRoundingMode(prevroundmode);
4786 
4787  if( c < noperands )
4788  {
4789  BMScopyMemoryArray(&resultants[c], &operands[c], noperands - c); /*lint !e776 !e866*/
4790  }
4791 
4792  return nreductions;
4793 }
void SCIPintervalSetEntire(SCIP_Real infinity, SCIP_INTERVAL *resultant)
void SCIPintervalReciprocal(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
SCIP_Bool SCIPintervalIsSubsetEQ(SCIP_Real infinity, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalUnify(SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
int SCIP_ROUNDMODE
Definition: intervalarith.h:64
void SCIPintervalCos(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
SCIP_Bool SCIPintervalIsEntire(SCIP_Real infinity, SCIP_INTERVAL operand)
#define SCIP_ROUND_NEAREST
SCIP_Real SCIPintervalPowerScalarIntegerSup(SCIP_Real operand1, int operand2)
#define infinity
Definition: gastrans.c:80
SCIP_Bool SCIPintervalIsPositiveInfinity(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalPowerScalarScalar(SCIP_INTERVAL *resultant, SCIP_Real operand1, SCIP_Real operand2)
void SCIPintervalMulSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSetRoundingMode(SCIP_ROUNDMODE roundmode)
void SCIPintervalQuad(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL xrng)
SCIP_Bool SCIPintervalIsNegativeInfinity(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalSub(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIP_Real SCIPintervalPowerScalarIntegerInf(SCIP_Real operand1, int operand2)
void SCIPintervalScalprodScalarsInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_Real *operand2)
#define FALSE
Definition: def.h:96
#define EPSISINT(x, eps)
Definition: def.h:223
SCIP_Real SCIPrelDiff(SCIP_Real val1, SCIP_Real val2)
Definition: misc.c:11096
void SCIPintervalAddScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalIntersectEps(SCIP_INTERVAL *resultant, SCIP_Real eps, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSolveUnivariateQuadExpression(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
#define TRUE
Definition: def.h:95
void SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real sqrcoeff, SCIP_Real lincoeff, SCIP_Real rhs, SCIP_INTERVAL xbnds)
void SCIPintervalPowerScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
static SCIP_Real negate(SCIP_Real x)
void SCIPintervalSet(SCIP_INTERVAL *resultant, SCIP_Real value)
#define EPSGE(x, y, eps)
Definition: def.h:215
#define SCIPdebugMessage
Definition: pub_message.h:96
void SCIPintervalAdd(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIP_ROUNDMODE SCIPintervalGetRoundingMode(void)
SCIP_Real SCIPnegateReal(SCIP_Real x)
Definition: misc.c:10269
void SCIPintervalAddSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalMulInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIP_VAR ** x
Definition: circlepacking.c:63
static SCIP_ROUNDMODE intervalGetRoundingMode(void)
real eps
void SCIPintervalScalprodScalars(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_Real *operand2)
SCIP_Bool SCIPintervalAreDisjoint(SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalDiv(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIP_Real inf
Definition: intervalarith.h:55
#define EPSN(x, eps)
Definition: def.h:218
void SCIPintervalAbs(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
#define SCIPerrorMessage
Definition: pub_message.h:64
interval arithmetics for provable bounds
#define SCIPdebugPrintf
Definition: pub_message.h:99
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalSolveBivariateQuadExpressionAllScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real ax, SCIP_Real ay, SCIP_Real axy, SCIP_Real bx, SCIP_Real by, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds, SCIP_INTERVAL ybnds)
void SCIPintervalScalprod(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_INTERVAL *operand2)
void SCIPintervalSetEmpty(SCIP_INTERVAL *resultant)
void SCIPintervalScalprodScalarsSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_Real *operand2)
SCIP_Real SCIPintervalNegateReal(SCIP_Real x)
void SCIPintervalSquareRoot(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
#define SCIP_ROUND_ZERO
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
internal miscellaneous methods
#define NULL
Definition: lpi_spx1.cpp:164
SCIP_Real SCIPintervalQuadUpperBound(SCIP_Real infinity, SCIP_Real a, SCIP_INTERVAL b_, SCIP_INTERVAL x)
SCIP_Real SCIPintervalGetSup(SCIP_INTERVAL interval)
void SCIPintervalMulScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalMulScalarSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
#define REALABS(x)
Definition: def.h:210
void SCIPintervalIntersect(SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
int SCIPintervalPropagateWeightedSum(SCIP_Real infinity, int noperands, SCIP_INTERVAL *operands, SCIP_Real *weights, SCIP_Real constant, SCIP_INTERVAL rhs, SCIP_INTERVAL *resultants, SCIP_Bool *infeasible)
#define SCIP_ROUND_UPWARDS
SCIP_Real sup
Definition: intervalarith.h:56
void SCIPintervalEntropy(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalLog(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalSolveUnivariateQuadExpressionPositive(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
void SCIPintervalSetRoundingModeDownwards(void)
void SCIPintervalSetRoundingModeToNearest(void)
static const double pi_d_u
#define SCIP_Bool
Definition: def.h:93
void SCIPintervalAddInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalPower(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
#define MAX(x, y)
Definition: tclique_def.h:92
void SCIPintervalPowerScalarInverse(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL basedomain, SCIP_Real exponent, SCIP_INTERVAL image)
SCIP_Real SCIPnextafter(SCIP_Real from, SCIP_Real to)
Definition: misc.c:9275
void SCIPintervalPowerScalarInteger(SCIP_INTERVAL *resultant, SCIP_Real operand1, int operand2)
#define EPSLE(x, y, eps)
Definition: def.h:213
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:136
SCIP_Bool SCIPintervalAreDisjointEps(SCIP_Real eps, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSubScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
#define SCIP_ROUND_DOWNWARDS
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 SCIPintervalSign(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalSetRoundingModeTowardsZero(void)
#define SCIP_REAL_MAX
Definition: def.h:187
void SCIPintervalSignPowerScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
#define EPSLT(x, y, eps)
Definition: def.h:212
#define SCIP_REAL_MIN
Definition: def.h:188
void SCIPintervalDivScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
#define EPSGT(x, y, eps)
Definition: def.h:214
SCIP_VAR ** b
Definition: circlepacking.c:65
static const double pi_d_l
void SCIPintervalSquare(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalMulScalarInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalMul(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSetRoundingModeUpwards(void)
public methods for message output
#define CALCB(y)
SCIP_VAR * a
Definition: circlepacking.c:66
SCIP_Bool SCIPintervalHasRoundingControl(void)
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
#define SCIP_Real
Definition: def.h:186
SCIP_VAR ** y
Definition: circlepacking.c:64
static void intervalSetRoundingMode(SCIP_ROUNDMODE roundmode)
void SCIPintervalAddVectors(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_INTERVAL *operand2)
#define CALCR(c, y)
common defines and data types used in all packages of SCIP
void SCIPintervalMax(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSolveUnivariateQuadExpressionNegative(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
void SCIPintervalMin(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSin(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalExp(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
#define EPSZ(x, eps)
Definition: def.h:216