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 */
89   	SCIP_Bool SCIPintervalHasRoundingControl(
90   	   void
91   	   )
92   	{
93   	   return TRUE;
94   	}
95   	
96   	/** sets rounding mode of floating point operations */
97   	static
98   	void intervalSetRoundingMode(
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
115  	SCIP_ROUNDMODE intervalGetRoundingMode(
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 */
141  	SCIP_Bool SCIPintervalHasRoundingControl(
142  	   void
143  	   )
144  	{
145  	   return TRUE;
146  	}
147  	
148  	/** sets rounding mode of floating point operations */
149  	static
150  	void intervalSetRoundingMode(
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
167  	SCIP_ROUNDMODE intervalGetRoundingMode(
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 */
193  	SCIP_Bool SCIPintervalHasRoundingControl(
194  	   void
195  	   )
196  	{
197  	   return TRUE;
198  	}
199  	
200  	/** sets rounding mode of floating point operations */
201  	static
202  	void intervalSetRoundingMode(
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
219  	SCIP_ROUNDMODE intervalGetRoundingMode(
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 */
239  	SCIP_Bool SCIPintervalHasRoundingControl(
240  	   void
241  	   )
242  	{
243  	   return FALSE;
244  	}
245  	
246  	/** sets rounding mode of floating point operations */ /*lint -e715*/
247  	static
248  	void intervalSetRoundingMode(
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
257  	SCIP_ROUNDMODE intervalGetRoundingMode(
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 */
268  	void SCIPintervalSetRoundingMode(
269  	   SCIP_ROUNDMODE        roundmode           /**< rounding mode to activate */
270  	   )
271  	{
272  	   intervalSetRoundingMode(roundmode);
273  	}
274  	
275  	/** gets current rounding mode of floating point operations */
276  	SCIP_ROUNDMODE SCIPintervalGetRoundingMode(
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
333  	SCIP_Real negate(
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 */
344  	void SCIPintervalSetRoundingModeDownwards(
345  	   void
346  	   )
347  	{
348  	   intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
349  	}
350  	
351  	/** sets rounding mode of floating point operations to upwards rounding */
352  	void SCIPintervalSetRoundingModeUpwards(
353  	   void
354  	   )
355  	{
356  	   intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
357  	}
358  	
359  	/** sets rounding mode of floating point operations to nearest rounding */
360  	void SCIPintervalSetRoundingModeToNearest(
361  	   void
362  	   )
363  	{
364  	   intervalSetRoundingMode(SCIP_ROUND_NEAREST);
365  	}
366  	
367  	/** sets rounding mode of floating point operations to towards zero rounding */
368  	void SCIPintervalSetRoundingModeTowardsZero(
369  	   void
370  	   )
371  	{
372  	   intervalSetRoundingMode(SCIP_ROUND_ZERO);
373  	}
374  	
375  	/** negates a number in a way that the compiler does not optimize it away */
376  	SCIP_Real SCIPintervalNegateReal(
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 */
405  	SCIP_Real SCIPintervalGetInf(
406  	   SCIP_INTERVAL         interval            /**< interval */
407  	   )
408  	{
409  	   return interval.inf;
410  	}
411  	
412  	/** returns supremum of interval */
413  	SCIP_Real SCIPintervalGetSup(
414  	   SCIP_INTERVAL         interval            /**< interval */
415  	   )
416  	{
417  	   return interval.sup;
418  	}
419  	
420  	/** stores given value as interval */
421  	void SCIPintervalSet(
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 */
433  	void SCIPintervalSetBounds(
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] */
447  	void SCIPintervalSetEmpty(
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 */
458  	SCIP_Bool SCIPintervalIsEmpty(
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] */
470  	void SCIPintervalSetEntire(
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 */
482  	SCIP_Bool SCIPintervalIsEntire(
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] */
491  	SCIP_Bool SCIPintervalIsPositiveInfinity(
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] */
500  	SCIP_Bool SCIPintervalIsNegativeInfinity(
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 */
509  	SCIP_Bool SCIPintervalIsSubsetEQ(
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 */
528  	SCIP_Bool SCIPintervalAreDisjoint(
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  	 */
541  	SCIP_Bool SCIPintervalAreDisjointEps(
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 */
557  	void SCIPintervalIntersect(
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  	 */
578  	void SCIPintervalIntersectEps(
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 */
609  	void SCIPintervalUnify(
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 */
636  	void SCIPintervalAddInf(
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  	{
643  	   assert(intervalGetRoundingMode() == SCIP_ROUND_DOWNWARDS);
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 */
663  	void SCIPintervalAddSup(
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  	{
670  	   assert(intervalGetRoundingMode() == SCIP_ROUND_UPWARDS);
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 */
690  	void SCIPintervalAdd(
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 */
706  	   intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
707  	   SCIPintervalAddInf(infinity, resultant, operand1, operand2);
708  	
709  	   /* compute supremum of result */
710  	   intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
711  	   SCIPintervalAddSup(infinity, resultant, operand1, operand2);
712  	
713  	   intervalSetRoundingMode(roundmode);
714  	}
715  	
716  	/** adds operand1 and scalar operand2 and stores result in resultant */
717  	void SCIPintervalAddScalar(
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  	   {
743  	      intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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  	   {
759  	      intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 */
767  	void SCIPintervalAddVectors(
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 */
781  	   intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
782  	   for( i = 0; i < length; ++i )
783  	   {
784  	      SCIPintervalAddInf(infinity, &resultant[i], operand1[i], operand2[i]);
785  	   }
786  	   /* compute supremums of result array */
787  	   intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 */
797  	void SCIPintervalSub(
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  	   {
823  	      intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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  	   {
837  	      intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 */
845  	void SCIPintervalSubScalar(
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 */
856  	void SCIPintervalMulInf(
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  	
867  	   assert(intervalGetRoundingMode() == SCIP_ROUND_DOWNWARDS);
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 */
916  	void SCIPintervalMulSup(
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  	
927  	   assert(intervalGetRoundingMode() == SCIP_ROUND_UPWARDS);
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 */
976  	void SCIPintervalMul(
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 */
992  	   intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
993  	   SCIPintervalMulInf(infinity, resultant, operand1, operand2);
994  	
995  	   /* compute supremum of result */
996  	   intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 */
1003 	void SCIPintervalMulScalarInf(
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 	{
1010 	   assert(intervalGetRoundingMode() == SCIP_ROUND_DOWNWARDS);
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 */
1059 	void SCIPintervalMulScalarSup(
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 	{
1066 	   assert(intervalGetRoundingMode() == SCIP_ROUND_UPWARDS);
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 */
1115 	void SCIPintervalMulScalar(
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 */
1143 	   intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
1144 	   SCIPintervalMulScalarInf(infinity, resultant, operand1, operand2);
1145 	
1146 	   /* compute supremum of result */
1147 	   intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
1148 	   SCIPintervalMulScalarSup(infinity, resultant, operand1, operand2);
1149 	
1150 	   intervalSetRoundingMode(roundmode);
1151 	}
1152 	
1153 	/** divides operand1 by operand2 and stores result in resultant */
1154 	void SCIPintervalDiv(
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 	   {
1190 	      intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 	   {
1199 	      intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 */
1208 	void SCIPintervalDivScalar(
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 	      {
1239 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 	      {
1252 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	      {
1267 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 	      {
1280 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 */
1311 	void SCIPintervalScalprod(
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 */
1329 	   intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 */
1339 	   intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 */
1351 	void SCIPintervalScalprodScalarsInf(
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 	
1362 	   assert(intervalGetRoundingMode() == SCIP_ROUND_DOWNWARDS);
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 */
1377 	void SCIPintervalScalprodScalarsSup(
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 	
1388 	   assert(intervalGetRoundingMode() == SCIP_ROUND_UPWARDS);
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 */
1403 	void SCIPintervalScalprodScalars(
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 */
1419 	   intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
1420 	   SCIPintervalScalprodScalarsInf(infinity, resultant, length, operand1, operand2);
1421 	   assert(resultant->sup == 0.0);
1422 	
1423 	   /* compute supremum of resultant */
1424 	   intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
1425 	   SCIPintervalScalprodScalarsSup(infinity, resultant, length, operand1, operand2);
1426 	
1427 	   intervalSetRoundingMode(roundmode);
1428 	}
1429 	
1430 	/** squares operand and stores result in resultant */
1431 	void SCIPintervalSquare(
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 	      {
1450 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
1451 	         resultant->inf = operand.sup * operand.sup;
1452 	      }
1453 	
1454 	      if( operand.inf <= -infinity )
1455 	         resultant->sup = infinity;
1456 	      else
1457 	      {
1458 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	      {
1468 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
1469 	         resultant->inf = operand.inf * operand.inf;
1470 	      }
1471 	
1472 	      if( operand.sup >= infinity )
1473 	         resultant->sup = infinity;
1474 	      else
1475 	      {
1476 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	
1490 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	 */
1503 	void SCIPintervalSquareRoot(
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 	 */
1564 	void SCIPintervalPower(
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 	 */
1610 	SCIP_Real SCIPintervalPowerScalarIntegerInf(
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();
1640 	      intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 	
1655 	      intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 	 */
1693 	SCIP_Real SCIPintervalPowerScalarIntegerSup(
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();
1723 	      intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	
1738 	      intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	 */
1768 	void SCIPintervalPowerScalarInteger(
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 	
1823 	      intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	 */
1860 	void SCIPintervalPowerScalarScalar(
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 	 */
1903 	void SCIPintervalPowerScalar(
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 	 */
2152 	void SCIPintervalPowerScalarInverse(
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 	 */
2233 	void SCIPintervalSignPowerScalar(
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 	      {
2303 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
2304 	         resultant->inf = operand1.inf * operand1.inf;
2305 	      }
2306 	      else
2307 	      {
2308 	         /* need upwards since we negate result of multiplication */
2309 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	      {
2323 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
2324 	         resultant->sup = operand1.sup * operand1.sup;
2325 	      }
2326 	      else
2327 	      {
2328 	         /* need downwards since we negate result of multiplication */
2329 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 	      {
2342 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
2343 	         resultant->inf =  SCIPnextafter(sqrt( operand1.inf), SCIP_REAL_MIN);
2344 	      }
2345 	      else
2346 	      {
2347 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	      {
2357 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
2358 	         resultant->sup =  SCIPnextafter(sqrt( operand1.sup), SCIP_REAL_MAX);
2359 	      }
2360 	      else
2361 	      {
2362 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	      {
2375 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
2376 	         resultant->inf =  SCIPnextafter(pow( operand1.inf, operand2), SCIP_REAL_MIN);
2377 	      }
2378 	      else
2379 	      {
2380 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	      {
2390 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
2391 	         resultant->sup =  SCIPnextafter(pow( operand1.sup, operand2), SCIP_REAL_MAX);
2392 	      }
2393 	      else
2394 	      {
2395 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	 */
2405 	void SCIPintervalReciprocal(
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 	      {
2431 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 	      {
2441 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	      {
2455 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
2456 	         resultant->inf = 1.0 / operand.sup;
2457 	      }
2458 	
2459 	      if( operand.inf <= -infinity )
2460 	         resultant->sup = infinity;
2461 	      else
2462 	      {
2463 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	 */
2478 	void SCIPintervalExp(
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 	
2512 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	
2534 	      assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	   {
2552 	      assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	 */
2562 	void SCIPintervalLog(
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 	
2591 	         assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	   {
2610 	      assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	   {
2624 	      assert(intervalGetRoundingMode() == SCIP_ROUND_NEAREST);
2625 	      resultant->sup = SCIPnextafter(log(operand.sup), SCIP_REAL_MAX);
2626 	   }
2627 	}
2628 	
2629 	/** stores minimum of operands in resultant */
2630 	void SCIPintervalMin(
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 */
2646 	void SCIPintervalMax(
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 */
2662 	void SCIPintervalAbs(
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 */
2700 	void SCIPintervalSin(
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 */
2742 	void SCIPintervalCos(
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 	
2771 	      assert(SCIPintervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 */
2790 	   SCIPintervalSetRoundingModeDownwards();
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( fmod(k, 2.0) != 0.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 */
2869 	void SCIPintervalSign(
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 */
2896 	void SCIPintervalEntropy(
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) */
2932 	   assert(SCIPintervalGetRoundingMode() == SCIP_ROUND_NEAREST);
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 	    */
2950 	   SCIPintervalSetRoundingModeUpwards();
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) */
2974 	   SCIPintervalSetRoundingModeToNearest();
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 	 */
2995 	SCIP_Real SCIPintervalQuadUpperBound(
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();
3027 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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();
3061 	      intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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$ */
3090 	void SCIPintervalQuad(
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 	 */
3125 	void SCIPintervalSolveUnivariateQuadExpressionPositive(
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 	 */
3172 	void SCIPintervalSolveUnivariateQuadExpressionNegative(
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 	 */
3208 	void SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(
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 	
3260 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	
3272 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 */
3300 	         intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
3301 	         delta = b*b + sqrcoeff*rhs;
3302 	         if( delta < 0.0 )
3303 	         {
3304 	            SCIPintervalSetEmpty(resultant);
3305 	         }
3306 	         else
3307 	         {
3308 	            intervalSetRoundingMode(SCIP_ROUND_NEAREST);
3309 	            z = SCIPnextafter(sqrt(delta), SCIP_REAL_MAX);
3310 	            intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	         {
3321 	            intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
3322 	            delta = b*b + sqrcoeff*rhs;
3323 	            intervalSetRoundingMode(SCIP_ROUND_NEAREST);
3324 	            z = SCIPnextafter(sqrt(delta), SCIP_REAL_MAX);
3325 	            intervalSetRoundingMode(SCIP_ROUND_UPWARDS);
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 	         {
3337 	            intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
3338 	            delta = b*b + sqrcoeff*rhs;
3339 	            intervalSetRoundingMode(SCIP_ROUND_NEAREST);
3340 	            z = SCIPnextafter(sqrt(delta), SCIP_REAL_MIN);
3341 	            intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 */
3352 	         intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
3353 	         delta = b*b + sqrcoeff * rhs;
3354 	         if( delta >= 0.0 )
3355 	         {
3356 	            /* let resultant = [0,-c/z] for now */
3357 	            intervalSetRoundingMode(SCIP_ROUND_NEAREST);
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 	             */
3362 	            intervalSetRoundingMode(SCIP_ROUND_DOWNWARDS);
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 	 */
3412 	void SCIPintervalSolveUnivariateQuadExpression(
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 	 */
3478 	void SCIPintervalQuadBivar(
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 	 */
3739 	void SCIPintervalSolveBivariateQuadExpressionAllScalar(
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 	 */
4625 	int SCIPintervalPropagateWeightedSum(
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();
4658 	   SCIPintervalSetRoundingModeDownwards();
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 	}
4794