35#define _USE_MATH_DEFINES
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
63#if defined(__GNUC__) && !defined(__INTEL_COMPILER)
65#pragma clang optimize off
67#pragma GCC optimize ("O0")
74#ifdef SCIP_ROUNDING_FE
83#define SCIP_ROUND_DOWNWARDS FE_DOWNWARD
84#define SCIP_ROUND_UPWARDS FE_UPWARD
85#define SCIP_ROUND_NEAREST FE_TONEAREST
86#define SCIP_ROUND_ZERO FE_TOWARDZERO
103 if( fesetround(roundmode) != 0 )
109 (void) fesetround(roundmode);
126#ifdef SCIP_ROUNDING_FP
135#define SCIP_ROUND_DOWNWARDS FP_RND_RM
136#define SCIP_ROUND_UPWARDS FP_RND_RP
137#define SCIP_ROUND_NEAREST FP_RND_RN
138#define SCIP_ROUND_ZERO FP_RND_RZ
155 if( write_rnd(roundmode) != 0 )
161 (void) write_rnd(roundmode);
178#ifdef SCIP_ROUNDING_MS
187#define SCIP_ROUND_DOWNWARDS RC_DOWN
188#define SCIP_ROUND_UPWARDS RC_UP
189#define SCIP_ROUND_NEAREST RC_NEAR
190#define SCIP_ROUND_ZERO RC_CHOP
207 if( (_controlfp(roundmode, _MCW_RC) & _MCW_RC) != roundmode )
213 (void) _controlfp(roundmode, _MCW_RC);
223 return _controlfp(0, 0) & _MCW_RC;
233#define SCIP_ROUND_DOWNWARDS 0
234#define SCIP_ROUND_UPWARDS 1
235#define SCIP_ROUND_NEAREST 2
236#define SCIP_ROUND_ZERO 3
252 SCIPerrorMessage(
"setting rounding mode not available - interval arithmetic is invalid!\n");
283#if defined(__GNUC__) && (defined(__i386__) || defined(__x86_64__))
297 __asm
volatile (
"fldl %1; fchs; fstpl %0" :
"=m" (
x) :
"m" (
x));
304#elif defined(_MSC_VER) && (defined(__INTEL_COMPILER) || !defined(_M_X64))
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
428 resultant->
inf = value;
429 resultant->
sup = value;
442 resultant->
inf = inf;
443 resultant->
sup = sup;
453 resultant->
inf = 1.0;
454 resultant->
sup = -1.0;
466 return operand.
sup < operand.
inf;
516 if( operand1.
inf > operand1.
sup )
520 if( operand2.
inf > operand2.
sup )
533 return (operand1.
sup < operand2.
inf) || (operand2.
sup < operand1.
inf);
547 if( operand1.
sup < operand2.
inf )
550 if( operand1.
inf > operand2.
sup )
588 if( operand1.
sup < operand2.
inf )
596 else if( operand1.
inf > operand2.
sup )
617 if( operand1.
inf > operand1.
sup )
620 *resultant = operand2;
624 if( operand2.
inf > operand2.
sup )
627 *resultant = operand1;
658 resultant->
inf = operand1.
inf + operand2.
inf;
685 resultant->
sup = operand1.
sup + operand2.
sup;
744 resultant->
inf = operand1.
inf + operand2;
760 resultant->
sup = operand1.
sup + operand2;
782 for(
i = 0;
i < length; ++
i )
788 for(
i = 0;
i < length; ++
i )
824 resultant->
inf = operand1.
inf - operand2.
sup;
838 resultant->
sup = operand1.
sup - operand2.
inf;
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));
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));
1017 if( operand1.
inf > 0 )
1019 else if( operand1.
inf < 0 )
1022 resultant->
inf = 0.0;
1027 if( operand1.
sup > 0 )
1029 else if( operand1.
sup < 0 )
1032 resultant->
inf = 0.0;
1034 else if( operand2 == 0.0 )
1036 resultant->
inf = 0.0;
1038 else if( operand2 > 0.0 )
1045 resultant->
inf = operand1.
inf * operand2;
1054 resultant->
inf = operand1.
sup * operand2;
1073 if( operand1.
sup > 0 )
1075 else if( operand1.
sup < 0 )
1078 resultant->
sup = 0.0;
1083 if( operand1.
inf > 0 )
1085 else if( operand1.
inf < 0 )
1088 resultant->
sup = 0.0;
1090 else if( operand2 == 0.0 )
1092 resultant->
sup = 0.0;
1094 else if( operand2 > 0.0 )
1101 resultant->
sup = operand1.
sup * operand2;
1110 resultant->
sup = operand1.
inf * operand2;
1127 if( operand2 == 1.0 )
1129 *resultant = operand1;
1133 if( operand2 == -1.0 )
1135 resultant->
inf = -operand1.
sup;
1136 resultant->
sup = -operand1.
inf;
1168 if( operand2.
inf <= 0.0 && operand2.
sup >= 0.0 )
1175 if( operand1.
inf == 0.0 && operand1.
sup == 0.0 )
1191 intmed.
inf = 1.0 / operand2.
sup;
1200 intmed.
sup = 1.0 / operand2.
inf;
1225 resultant->
inf = 0.0;
1226 resultant->
sup = 0.0;
1228 else if( operand2 > 0.0 )
1240 resultant->
inf = operand1.
inf / operand2;
1253 resultant->
sup = operand1.
sup / operand2;
1256 else if( operand2 < 0.0 )
1268 resultant->
inf = operand1.
sup / operand2;
1281 resultant->
sup = operand1.
inf / operand2;
1286 if( operand1.
inf >= 0 )
1292 else if( operand1.
sup <= 0 )
1325 resultant->
inf = 0.0;
1326 resultant->
sup = 0.0;
1364 resultant->
inf = 0.0;
1390 resultant->
sup = 0.0;
1415 resultant->
inf = 0.0;
1416 resultant->
sup = 0.0;
1444 if( operand.
sup <= 0.0 )
1451 resultant->
inf = operand.
sup * operand.
sup;
1459 resultant->
sup = operand.
inf * operand.
inf;
1462 else if( operand.
inf >= 0.0 )
1469 resultant->
inf = operand.
inf * operand.
inf;
1477 resultant->
sup = operand.
sup * operand.
sup;
1482 resultant->
inf = 0.0;
1491 x = operand.
inf * operand.
inf;
1492 y = operand.
sup * operand.
sup;
1512 if( operand.
sup < 0.0 )
1518 if( operand.
inf == operand.
sup )
1530 tmp = sqrt(operand.
inf);
1538 if( operand.
inf <= 0.0 )
1539 resultant->
inf = 0.0;
1575 if( operand2.
inf == operand2.
sup )
1584 if( operand1.
sup == 0.0 )
1586 if( operand2.
inf <= 0.0 && operand2.
sup >= 0.0 )
1620 if( operand1 == 0.0 )
1630 if( operand1 == 1.0 || operand2 == 0 )
1652 n = (
unsigned int)operand2;
1703 if( operand1 == 0.0 )
1713 if( operand1 == 1.0 || operand2 == 0 )
1735 n = (
unsigned int)operand2;
1778 if( operand1 == 0.0 )
1794 if( operand1 == 1.0 || operand2 == 0 )
1812 SCIP_Real result_sup;
1813 SCIP_Real result_inf;
1821 n = (
unsigned int)operand2;
1834 result_sup = result_sup * z_sup;
1842 z_sup = z_sup * z_sup;
1848 resultant->
inf = result_inf;
1849 resultant->
sup = result_sup;
1870 if( operand1 == 0.0 )
1886 if( operand1 == 1.0 || operand2 == 0 )
1893 result = pow(operand1, operand2);
1921 if( operand1.
inf < 0.0 )
1924 resultant->
inf = 0.0;
1925 if( operand1.
sup > 0.0 )
1928 resultant->
sup = 0.0;
1932 if( operand2 == 0.0 )
1934 if( operand1.
inf == 0.0 && operand1.
sup == 0.0 )
1936 resultant->
inf = 0.0;
1937 resultant->
sup = 0.0;
1939 else if( operand1.
inf <= 0.0 || operand1.
sup >= 0.0 )
1941 resultant->
inf = 0.0;
1942 resultant->
sup = 1.0;
1946 resultant->
inf = 1.0;
1947 resultant->
sup = 1.0;
1952 if( operand2 == 1.0 )
1955 *resultant = operand1;
1959 op2isint = (ceil(operand2) == operand2);
1961 if( !op2isint && operand1.
inf < 0.0 )
1964 if( operand1.
sup < operand1.
inf )
1971 if( operand1.
inf >= 0.0 )
1973 if( operand2 >= 0.0 )
1978 else if( operand1.
inf > 0.0 )
1984 resultant->
inf = 0.0;
1988 else if( operand1.
sup > 0.0 )
1994 resultant->
sup = 0.0;
1999 resultant->
inf = 0.0;
2000 else if( operand1.
sup == 0.0 )
2004 if( ceil(operand2/2) == operand2/2 )
2016 if( operand1.
inf == 0.0 )
2025 else if( operand1.
sup <= 0.0 )
2028 if( operand2 >= 0.0 && ceil(operand2/2) == operand2/2 )
2042 else if( operand2 <= 0.0 && ceil(operand2/2) != operand2/2 )
2047 resultant->
inf = 0.0;
2048 else if( operand1.
sup == 0.0 )
2056 resultant->
sup = 0.0;
2057 else if( operand1.
inf == 0.0 )
2063 else if( operand2 >= 0.0 )
2080 resultant->
inf = 0.0;
2081 else if( operand1.
inf == 0.0 )
2088 resultant->
sup = 0.0;
2089 else if( operand1.
sup == 0.0 )
2100 if( operand2 >= 0.0 && operand2/2 == ceil(operand2/2) )
2103 resultant->
inf = 0.0;
2109 else if( operand2 <= 0.0 && ceil(operand2/2) == operand2/2 )
2114 resultant->
inf = 0.0;
2118 else if( operand2 >= 0.0 )
2167 if( exponent == 0.0 )
2170 if( image.
inf <= 1.0 && image.
sup >= 1.0 )
2173 *resultant = basedomain;
2175 else if( image.
inf <= 0.0 && image.
sup >= 0.0 )
2199 if( image.
sup >= 0.0 )
2203 if( basedomain.
inf <= -resultant->
inf &&
EPSISINT(exponent, 0.0) && (
int)exponent % 2 == 0 )
2205 if( basedomain.
sup < resultant->
inf )
2217 if( image.
inf < 0.0 && basedomain.
inf < 0.0 &&
EPSISINT(exponent, 0.0) && ((
int)exponent % 2 != 0) )
2252 if( operand1.
inf < 0.0 )
2255 resultant->
inf = 0.0;
2256 if( operand1.
sup > 0.0 )
2259 resultant->
sup = 0.0;
2263 if( operand2 == 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;
2271 resultant->
inf = 1.0;
2273 if( operand1.
sup < 0.0 )
2274 resultant->
sup = -1.0;
2275 else if( operand1.
sup == 0.0 )
2276 resultant->
sup = 0.0;
2278 resultant->
sup = 1.0;
2283 if( operand2 == 1.0 )
2285 *resultant = operand1;
2291 if( operand2 == 2.0 )
2301 else if( operand1.
inf > 0.0 )
2304 resultant->
inf = operand1.
inf * operand1.
inf;
2321 else if( operand1.
sup > 0.0 )
2324 resultant->
sup = operand1.
sup * operand1.
sup;
2334 else if( operand2 == 0.5 )
2340 else if( operand1.
inf >= 0.0 )
2355 else if( operand1.
sup > 0.0 )
2373 else if( operand1.
inf > 0.0 )
2388 else if( operand1.
sup > 0.0 )
2416 if( operand.
inf == 0.0 && operand.
sup == 0.0 )
2425 if( operand.
inf >= 0.0 )
2428 resultant->
inf = 0.0;
2432 resultant->
inf = 1.0 / operand.
sup;
2436 resultant->
sup = 0.0;
2437 else if( operand.
inf == 0.0 )
2442 resultant->
sup = 1.0 / operand.
inf;
2447 else if( operand.
sup <= 0.0 )
2450 resultant->
inf = 0.0;
2451 else if( operand.
sup == 0.0 )
2456 resultant->
inf = 1.0 / operand.
sup;
2464 resultant->
sup = 1.0 / operand.
inf;
2489 resultant->
inf = 0.0;
2490 resultant->
sup = 0.0;
2501 if( operand.
inf == operand.
sup )
2503 if( operand.
inf == 0.0 )
2505 resultant->
inf = 1.0;
2506 resultant->
sup = 1.0;
2513 tmp = exp(operand.
inf);
2524 resultant->
inf = 0.0;
2526 else if( operand.
inf == 0.0 )
2528 resultant->
inf = 1.0;
2535 tmp = exp(operand.
inf);
2546 else if( operand.
sup == 0.0 )
2548 resultant->
sup = 1.0;
2574 if( operand.
sup <= 0.0 )
2580 if( operand.
inf == operand.
sup )
2582 if( operand.
sup == 1.0 )
2584 resultant->
inf = 0.0;
2585 resultant->
sup = 0.0;
2592 tmp = log(operand.
inf);
2600 if( operand.
inf <= 0.0 )
2604 else if( operand.
inf == 1.0 )
2606 resultant->
inf = 0.0;
2618 else if( operand.
sup == 1.0 )
2620 resultant->
sup = 0.0;
2671 if( operand.
inf <= 0.0 && operand.
sup >= 0.0)
2673 resultant->
inf = 0.0;
2676 else if( operand.
inf > 0.0 )
2678 *resultant = operand;
2682 resultant->
inf = -operand.
sup;
2683 resultant->
sup = -operand.
inf;
2692static const double pi_d_l = (3373259426.0 + 273688.0 / (1<<21)) / (1<<30);
2693static const double pi_d_u = (3373259426.0 + 273689.0 / (1<<21)) / (1<<30);
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))
2724 tmp = -resultant->
sup;
2725 resultant->
sup = -resultant->
inf;
2726 resultant->
inf = tmp;
2736 resultant->
inf = 0.0;
2738 resultant->
sup = 0.0;
2767 if( operand.
inf == operand.
sup )
2772 tmp = cos(operand.
inf);
2781 if( operand.
sup > 1e12 || operand.
inf < -1e12 )
2791 negwidth = operand.
inf - operand.
sup;
2792 if( -negwidth >= 2.0*
pi_d_l )
2822 resultant->
inf =
MAX(-1.0, resultant->
inf);
2823 if( operand.
inf == 0.0 )
2824 resultant->
sup = 1.0;
2828 resultant->
sup =
MIN( 1.0, resultant->
sup);
2835 resultant->
inf = -1.0;
2836 if( operand.
inf == 0.0 )
2837 resultant->
sup = 1.0;
2843 cinf = cos(operand.
inf);
2844 csup = cos(operand.
sup);
2846 resultant->
sup =
MIN(1.0, resultant->
sup);
2856 if( fmod(k, 2.0) != 0.0 )
2858 SCIP_Real tmp = -resultant->
sup;
2859 resultant->
sup = -resultant->
inf;
2860 resultant->
inf = tmp;
2878 if( operand.
sup < 0.0 )
2880 resultant->
inf = -1.0;
2881 resultant->
sup = -1.0;
2883 else if( operand.
inf >= 0.0 )
2885 resultant->
inf = 1.0;
2886 resultant->
sup = 1.0;
2890 resultant->
inf = -1.0;
2891 resultant->
sup = 1.0;
2904 SCIP_Real infcand1 = 0.0;
2905 SCIP_Real infcand2 = 0.0;
2906 SCIP_Real supcand1 = 0.0;
2907 SCIP_Real supcand2 = 0.0;
2916 if( operand.
sup < 0.0 )
2923 if( operand.
sup == 0.0 )
2933 if( operand.
inf > 0.0 )
2935 loginf = log(operand.
inf);
2942 logsup = log(operand.
sup);
2951 if( operand.
inf > 0.0 )
2976 inf =
MIN(infcand1, infcand2);
2979 if( operand.
inf <= extr && extr <= operand.
sup )
2982 sup =
MAX3(supcand1, supcand2, extr);
2985 sup =
MAX(supcand1, supcand2);
3024 SCIP_Real cand1, cand2, cand3, cand4;
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));
3063 u =
MAX(
x.inf * (
a*
x.inf +
b),
x.sup * (
a*
x.sup +
b));
3066 if( t >
x.inf &&
negate(2*
a)*
x.sup >
b && s*t > u )
3083 return MAX(cand1, cand2);
3105 if( sqrcoeff == 0.0 )
3114 lincoeff.
inf = -lincoeff.
sup;
3115 lincoeff.
sup = -tmp;
3139 resultant->
inf = 0.0;
3188 lincoeff.
inf = -lincoeff.
sup;
3189 lincoeff.
sup = -tmp;
3197 tmp = resultant->
inf;
3198 resultant->
inf = -resultant->
sup;
3199 resultant->
sup = -tmp;
3226 if( sqrcoeff == 0.0 )
3239 if( lincoeff <= 0.0 && rhs > 0.0 )
3245 if( lincoeff >= 0.0 && rhs <= 0.0 )
3249 resultant->
sup = xbnds.
sup;
3255 if( lincoeff < 0.0 && rhs <= 0.0 )
3261 resultant->
sup = rhs / lincoeff;
3262 if( xbnds.
sup < resultant->
sup )
3263 resultant->
sup = xbnds.
sup;
3273 resultant->
inf = rhs / lincoeff;
3274 if( resultant->
inf < xbnds.
inf )
3275 resultant->
inf = xbnds.
inf;
3277 resultant->
sup = xbnds.
sup;
3285 resultant->
inf = 0.0;
3296 if( lincoeff >= 0.0 )
3301 delta =
b*
b + sqrcoeff*rhs;
3313 if( sqrcoeff < 0.0 )
3319 if( sqrcoeff < 0.0 )
3322 delta =
b*
b + sqrcoeff*rhs;
3335 if( sqrcoeff > 0.0 )
3338 delta =
b*
b + sqrcoeff*rhs;
3343 resultant->
inf = z / sqrcoeff;
3353 delta =
b*
b + sqrcoeff * rhs;
3367 if( sqrcoeff > 0.0 )
3376 if( xbnds.
sup < zdiva )
3380 else if( xbnds.
inf > resultant->
sup )
3383 resultant->
inf = zdiva;
3433 if( sqrcoeff.
inf == 0.0 && sqrcoeff.
sup == 0.0 && (lincoeff.
inf > 0.0 || lincoeff.
sup < 0.0) )
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);
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);
3444 if( xbnds.
sup >= 0 )
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);
3456 if( xbnds.
inf <= 0.0 )
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);
3524 denom = 4.0 * ax * ay - axy * axy;
3527 x = (axy * by - 2.0 * ay * bx) / denom;
3528 y = (axy * bx - 2.0 * ax * by) / denom;
3531 val = (axy * bx * by - ay * bx * bx - ax * by * by) / denom;
3532 minval =
MIN(val, minval);
3533 maxval =
MAX(val, maxval);
3536 else if(
REALABS(2.0 * ay * bx - axy * by) <= 1e-9 )
3544 val = -ay * bx * bx / (axy * axy);
3545 minval =
MIN(val, minval);
3546 maxval =
MAX(val, maxval);
3559 else if( ax == 0.0 )
3565 else if( bx + axy * ybnds.
inf < 0.0 )
3569 minval =
MIN(val, minval);
3570 maxval =
MAX(val, maxval);
3574 else if( bx + axy * ybnds.
sup < 0.0 )
3578 minval =
MIN(val, minval);
3579 maxval =
MAX(val, maxval);
3591 minval =
MIN(tmp.
inf, minval);
3592 maxval =
MAX(tmp.
sup, maxval);
3602 else if( ax == 0.0 )
3608 else if( bx + axy * ybnds.
inf > 0.0 )
3612 minval =
MIN(val, minval);
3613 maxval =
MAX(val, maxval);
3617 else if( bx + axy * ybnds.
sup > 0.0 )
3621 minval =
MIN(val, minval);
3622 maxval =
MAX(val, maxval);
3634 minval =
MIN(tmp.
inf, minval);
3635 maxval =
MAX(tmp.
sup, maxval);
3645 else if( ay == 0.0 )
3651 else if( by + axy * xbnds.
inf < 0.0 )
3655 minval =
MIN(val, minval);
3656 maxval =
MAX(val, maxval);
3660 else if( by + axy * xbnds.
sup < 0.0 )
3664 minval =
MIN(val, minval);
3665 maxval =
MAX(val, maxval);
3677 minval =
MIN(tmp.
inf, minval);
3678 maxval =
MAX(tmp.
sup, maxval);
3688 else if( ay == 0.0 )
3694 else if( by + axy * xbnds.
inf > 0.0 )
3698 minval =
MIN(val, minval);
3699 maxval =
MAX(val, maxval);
3703 else if( by + axy * xbnds.
sup > 0.0 )
3707 minval =
MIN(val, minval);
3708 maxval =
MAX(val, maxval);
3720 minval =
MIN(tmp.
inf, minval);
3721 maxval =
MAX(tmp.
sup, maxval);
3724 minval -= 1e-10 *
REALABS(minval);
3725 maxval += 1e-10 *
REALABS(maxval);
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);
3773 if( xbnds.
sup >= 0.0 )
3782 if( xbnds.
inf < 0.0 )
3848 SCIP_Real minvalleft;
3849 SCIP_Real maxvalleft;
3850 SCIP_Real minvalright;
3851 SCIP_Real maxvalright;
3855 SCIP_Real rcoef_const;
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);
3868#define CALCB(y) ((bx + axy * (y)) / (2.0 * sqrtax))
3869#define CALCR(c,y) (rcoef_const + (c) + (rcoef_y + rcoef_yy * (y)) * (y))
3880 if(
EPSN(ub, 1e-9) )
3920 if( !
EPSZ(ay, 1e-9) && axy * axy >= 4.0 * ax * ay )
3941 else if( !
EPSZ(ay, 1e-9) )
3948 minvalleft = -by / 2.0;
3949 maxvalleft = -by / 2.0;
3970 minvalleft =
MIN(-sqrtc -
b, minvalleft);
3971 maxvalright =
MAX( sqrtc -
b, maxvalright);
3984 maxvalleft =
MAX(-sqrtc -
b, maxvalleft);
3985 minvalright =
MIN( sqrtc -
b, minvalright);
3994 if( !
EPSZ(ay, 1e-9) && axy * axy >= 4.0 * ax * ay )
4015 else if( !
EPSZ(ay, 1e-9) )
4024 minvalright =
MIN(minvalright, -by / 2.0);
4025 maxvalright =
MAX(maxvalright, -by / 2.0);
4044 minvalleft =
MIN(-sqrtc -
b, minvalleft);
4045 maxvalright =
MAX( sqrtc -
b, maxvalright);
4058 maxvalleft =
MAX(-sqrtc -
b, maxvalleft);
4059 minvalright =
MIN( sqrtc -
b, minvalright);
4066 if( !
EPSZ(ay, 1e-9) )
4068 if(
REALABS(axy*axy - 4.0*ax*ay) > 1e-9 )
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) )
4077 sqrtterm = sqrt(
MAX(sqrtterm, 0.0));
4079 ymin = axy * ay * bx - 2.0 * ax * ay * by - sqrtterm;
4081 ymin /= 4.0 * ax * ay - axy * axy;
4083 if( ymin > ybnds.
inf && ymin < ybnds.
sup )
4096 minvalleft =
MIN(-sqrtc -
b, minvalleft);
4097 maxvalright =
MAX( sqrtc -
b, maxvalright);
4102 ymin = axy * ay * bx - 2.0 * ax * ay * by + sqrtterm;
4104 ymin /= 4.0 * ax * ay - axy * axy;
4106 if( ymin > ybnds.
inf && ymin < ybnds.
sup )
4119 minvalleft =
MIN(-sqrtc -
b, minvalleft);
4120 maxvalright =
MAX( sqrtc -
b, maxvalright);
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) )
4131 sqrtterm = sqrt(
MAX(sqrtterm, 0.0));
4133 ymin = axy * ay * bx - 2.0 * ax * ay * by - sqrtterm;
4135 ymin /= 4.0 * ax * ay - axy * axy;
4137 if( ymin > ybnds.
inf && ymin < ybnds.
sup )
4150 maxvalleft =
MAX(-sqrtc -
b, maxvalleft);
4151 minvalright =
MIN( sqrtc -
b, minvalright);
4156 ymin = axy * ay * bx - 2.0 * ax * ay * by + sqrtterm;
4158 ymin /= 4.0 * ax * ay - axy * axy;
4160 if( ymin > ybnds.
inf && ymin < ybnds.
sup )
4173 maxvalleft =
MAX(-sqrtc -
b, maxvalleft);
4174 minvalright =
MIN( sqrtc -
b, minvalright);
4180 else if(
REALABS(2.0 * ay * bx - axy * by) > 1e-9 )
4184 ymin = - (4.0 * ay * bx * by - axy * by * by + 4.0 * axy * ay * rhs.
sup);
4186 ymin /= 2.0 * ay * bx - axy * by;
4188 if( ymin > ybnds.
inf && ymin < ybnds.
sup )
4201 minvalleft =
MIN(-sqrtc -
b, minvalleft);
4202 maxvalright =
MAX( sqrtc -
b, maxvalright);
4209 ymin = - (4.0 * ay * bx * by - axy * by * by + 4.0 * axy * ay * rhs.
inf);
4211 ymin /= 2.0 * ay * bx - axy * by;
4213 if( ymin > ybnds.
inf && ymin < ybnds.
sup )
4226 maxvalleft =
MAX(-sqrtc -
b, maxvalleft);
4227 minvalright =
MIN( sqrtc -
b, minvalright);
4258 if( ybnds.
sup >= 0.0 )
4267 minvalleft =
MIN(minvalleft, -
b);
4268 maxvalleft =
MAX(maxvalleft, -
b);
4269 minvalright =
MIN(minvalright, -
b);
4270 maxvalright =
MAX(maxvalright, -
b);
4275 minvalleft =
MIN(minvalleft, -
b);
4276 maxvalleft =
MAX(maxvalleft, -
b);
4277 minvalright =
MIN(minvalright, -
b);
4278 maxvalright =
MAX(maxvalright, -
b);
4301 if( ybnds.
inf < 0.0 )
4311 minvalleft =
MIN(minvalleft, -
b);
4312 maxvalleft =
MAX(maxvalleft, -
b);
4313 minvalright =
MIN(minvalright, -
b);
4314 maxvalright =
MAX(maxvalright, -
b);
4333 minvalleft =
MIN(minvalleft, -
b);
4334 maxvalleft =
MAX(maxvalleft, -
b);
4335 minvalright =
MIN(minvalright, -
b);
4336 maxvalright =
MAX(maxvalright, -
b);
4400 if(
EPSGE(-bx / axy, ybnds.
inf, 1e-9) &&
EPSLE(-bx / axy, ybnds.
sup, 1e-9) )
4410 if( xbnds.
inf < 0.0 && xbnds.
sup > 0.0 )
4428 if( lincoef.
inf == 0.0 && lincoef.
sup == 0.0 )
4431 if( myrhs.
inf <= 0.0 && myrhs.
sup >= 0.0 )
4436 else if( xbnds.
inf >= 0.0 )
4467 if( bx + axy * (axy > 0.0 ? ybnds.
inf : ybnds.
sup) > 0.0 )
4477 if(
EPSZ(ay, 1e-9) )
4479 else if( ay * axy < 0.0 )
4484 val = (
c - ay * ybnds.
inf * ybnds.
inf - by * ybnds.
inf) / (bx + axy * ybnds.
inf);
4485 minval =
MIN(val, minval);
4491 if(
EPSZ(ay, 1e-9) )
4492 minval =
MIN(minval, -by / axy);
4493 else if( ay * axy > 0.0 )
4498 val = (
c - ay * ybnds.
sup * ybnds.
sup - by * ybnds.
sup) / (bx + axy * ybnds.
sup);
4499 minval =
MIN(val, minval);
4502 if( !
EPSZ(ay, 1e-9) )
4504 d = ay * (ay * bx * bx - axy * (bx * by + axy *
c));
4505 if( !
EPSN(d, 1e-9) )
4507 ymin = -ay * bx + sqrt(
MAX(d, 0.0));
4510 if( ymin > ybnds.
inf && ymin < ybnds.
sup )
4512 assert(bx + axy * ymin != 0.0);
4514 val = (
c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4515 minval =
MIN(val, minval);
4518 ymin = -ay * bx - sqrt(
MAX(d, 0.0));
4521 if(ymin > ybnds.
inf && ymin < ybnds.
sup )
4523 assert(bx + axy * ymin != 0.0);
4525 val = (
c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4526 minval =
MIN(val, minval);
4537 if( bx + axy * (axy > 0.0 ? ybnds.
inf : ybnds.
sup) > 0.0 )
4547 if(
EPSZ(ay, 1e-9) )
4549 else if( ay * axy > 0.0 )
4554 val = (
c - ay * ybnds.
inf * ybnds.
inf - by * ybnds.
inf) / (bx + axy * ybnds.
inf);
4555 maxval =
MAX(val, maxval);
4561 if(
EPSZ(ay, 1e-9) )
4562 maxval =
MAX(maxval, -by / axy);
4563 else if( ay * axy < 0.0 )
4568 val = (
c - ay * ybnds.
sup * ybnds.
sup - by * ybnds.
sup) / (bx + axy * ybnds.
sup);
4569 maxval =
MAX(val, maxval);
4572 if( !
EPSZ(ay, 1e-9) )
4574 d = ay * (ay * bx * bx - axy * (bx * by + axy *
c));
4575 if( !
EPSN(d, 1e-9) )
4577 ymin = ay * bx + sqrt(
MAX(d, 0.0));
4580 if( ymin > ybnds.
inf && ymin < ybnds.
sup )
4582 assert(bx + axy * ymin != 0.0);
4583 val = (
c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4584 maxval =
MAX(val, maxval);
4587 ymin = ay * bx - sqrt(
MAX(d, 0.0));
4590 if( ymin > ybnds.
inf && ymin < ybnds.
sup )
4592 assert(bx + axy * ymin != 0.0);
4593 val = (
c - ay * ymin * ymin - by * ymin) / (bx + axy * ymin);
4594 maxval =
MAX(val, maxval);
4605 resultant->
inf = minval - 1e-10 *
REALABS(minval);
4609 resultant->
sup = maxval + 1e-10 *
REALABS(maxval);
4633 SCIP_Bool* infeasible
4638 SCIP_Real minlinactivity;
4639 SCIP_Real maxlinactivity;
4640 int minlinactivityinf;
4641 int maxlinactivityinf;
4642 int nreductions = 0;
4651 *infeasible =
FALSE;
4660 minlinactivity = constant;
4661 maxlinactivity = -constant;
4662 minlinactivityinf = 0;
4663 maxlinactivityinf = 0;
4665 SCIPdebugMessage(
"reverse prop with %d children: %.20g", noperands, constant);
4670 for(
c = 0;
c < noperands; ++
c )
4672 childbounds = operands[
c];
4685 ++maxlinactivityinf;
4689 maxlinactivity -= resultants[
c].
sup;
4693 ++minlinactivityinf;
4697 minlinactivity += resultants[
c].
inf;
4700 maxlinactivity = -maxlinactivity;
4703 minlinactivityinf ? -
infinity : minlinactivity,
4704 maxlinactivityinf ?
infinity : maxlinactivity,
4708 if( (minlinactivityinf >= 2 || rhs.
sup >=
infinity) && (maxlinactivityinf >= 2 || rhs.
inf <= -
infinity) )
4714 for(
c = 0;
c < noperands; ++
c )
4724 if( resultants[
c].inf <= -
infinity && minlinactivityinf <= 1 )
4726 assert(minlinactivityinf == 1);
4729 else if( minlinactivityinf == 0 )
4741 if( resultants[
c].sup >=
infinity && maxlinactivityinf <= 1 )
4743 assert(maxlinactivityinf == 1);
4744 childbounds.
inf = rhs.
inf - maxlinactivity;
4746 else if( maxlinactivityinf == 0 )
4748 childbounds.
inf = rhs.
inf - maxlinactivity + resultants[
c].
sup;
4780 if( resultants[
c].inf != operands[
c].inf || resultants[
c].sup != operands[
c].sup )
common defines and data types used in all packages of SCIP
SCIP_Real SCIPnextafter(SCIP_Real from, SCIP_Real to)
SCIP_Real SCIPrelDiff(SCIP_Real val1, SCIP_Real val2)
void SCIPintervalMulSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSetRoundingModeUpwards(void)
void SCIPintervalIntersectEps(SCIP_INTERVAL *resultant, SCIP_Real eps, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalAddSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalMulScalarInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalSetRoundingModeDownwards(void)
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
SCIP_Real SCIPintervalQuadUpperBound(SCIP_Real infinity, SCIP_Real a, SCIP_INTERVAL b_, SCIP_INTERVAL x)
SCIP_Bool SCIPintervalIsEntire(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalScalprodScalarsInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_Real *operand2)
void SCIPintervalSub(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSetEntire(SCIP_Real infinity, SCIP_INTERVAL *resultant)
SCIP_Bool SCIPintervalIsPositiveInfinity(SCIP_Real infinity, SCIP_INTERVAL operand)
SCIP_Real SCIPintervalPowerScalarIntegerSup(SCIP_Real operand1, int operand2)
void SCIPintervalSquareRoot(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalMulInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalScalprodScalarsSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_Real *operand2)
SCIP_ROUNDMODE SCIPintervalGetRoundingMode(void)
void SCIPintervalSetRoundingModeToNearest(void)
void SCIPintervalPower(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSolveUnivariateQuadExpression(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
void SCIPintervalSignPowerScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_Real SCIPintervalPowerScalarIntegerInf(SCIP_Real operand1, int operand2)
void SCIPintervalSetRoundingMode(SCIP_ROUNDMODE roundmode)
void SCIPintervalScalprodScalars(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_Real *operand2)
void SCIPintervalUnify(SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalAbs(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalSubScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalSetRoundingModeTowardsZero(void)
void SCIPintervalCos(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalSolveBivariateQuadExpressionAllScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real ax, SCIP_Real ay, SCIP_Real axy, SCIP_Real bx, SCIP_Real by, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds, SCIP_INTERVAL ybnds)
void SCIPintervalMax(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalIntersect(SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSquare(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalSet(SCIP_INTERVAL *resultant, SCIP_Real value)
SCIP_Bool SCIPintervalIsSubsetEQ(SCIP_Real infinity, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalPowerScalarInverse(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL basedomain, SCIP_Real exponent, SCIP_INTERVAL image)
SCIP_Bool SCIPintervalHasRoundingControl(void)
void SCIPintervalSin(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalSolveUnivariateQuadExpressionPositive(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
void SCIPintervalAddInf(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalMin(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalAddVectors(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_INTERVAL *operand2)
SCIP_Bool SCIPintervalIsNegativeInfinity(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalMulScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalReciprocal(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalPowerScalarInteger(SCIP_INTERVAL *resultant, SCIP_Real operand1, int operand2)
void SCIPintervalEntropy(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalMul(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalDiv(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalPowerScalarScalar(SCIP_INTERVAL *resultant, SCIP_Real operand1, SCIP_Real operand2)
void SCIPintervalQuad(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL xrng)
void SCIPintervalSign(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalPowerScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalLog(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalAddScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalMulScalarSup(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalDivScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_Bool SCIPintervalAreDisjointEps(SCIP_Real eps, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSolveUnivariateQuadExpressionNegative(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
SCIP_Real SCIPintervalGetSup(SCIP_INTERVAL interval)
void SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real sqrcoeff, SCIP_Real lincoeff, SCIP_Real rhs, SCIP_INTERVAL xbnds)
void SCIPintervalQuadBivar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real ax, SCIP_Real ay, SCIP_Real axy, SCIP_Real bx, SCIP_Real by, SCIP_INTERVAL xbnds, SCIP_INTERVAL ybnds)
SCIP_Real SCIPintervalNegateReal(SCIP_Real x)
SCIP_Bool SCIPintervalAreDisjoint(SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
int SCIPintervalPropagateWeightedSum(SCIP_Real infinity, int noperands, SCIP_INTERVAL *operands, SCIP_Real *weights, SCIP_Real constant, SCIP_INTERVAL rhs, SCIP_INTERVAL *resultants, SCIP_Bool *infeasible)
void SCIPintervalExp(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalAdd(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalScalprod(SCIP_Real infinity, SCIP_INTERVAL *resultant, int length, SCIP_INTERVAL *operand1, SCIP_INTERVAL *operand2)
void SCIPintervalSetEmpty(SCIP_INTERVAL *resultant)
assert(minobj< SCIPgetCutoffbound(scip))
static const double pi_d_l
#define SCIP_ROUND_NEAREST
static SCIP_Real negate(SCIP_Real x)
static SCIP_ROUNDMODE intervalGetRoundingMode(void)
#define SCIP_ROUND_UPWARDS
static const double pi_d_u
#define SCIP_ROUND_DOWNWARDS
static void intervalSetRoundingMode(SCIP_ROUNDMODE roundmode)
interval arithmetics for provable bounds
#define BMScopyMemoryArray(ptr, source, num)
SCIP_Real SCIPnegateReal(SCIP_Real x)
internal miscellaneous methods
public methods for message output