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