Scippy

SCIP

Solving Constraint Integer Programs

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