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