1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file intervalarith.c
17 * @brief unit tests to check that bug1861 doesn't appear
18 * @author Stefan Vigerske
19 */
20
21 #include "scip/intervalarith.h"
22 #include "include/scip_test.h"
23 #include <stdio.h>
24
Test(intervalarith,issue1861)25 Test(intervalarith, issue1861)
26 {
27 SCIP_Real infinity;
28 SCIP_INTERVAL resultant;
29 SCIP_Real ax;
30 SCIP_Real ay;
31 SCIP_Real axy;
32 SCIP_Real bx;
33 SCIP_Real by;
34 SCIP_INTERVAL rhs;
35 SCIP_INTERVAL xbnds;
36 SCIP_INTERVAL ybnds;
37
38 infinity = 1e43;
39 ax = 12;
40 ay = 27;
41 axy = -36;
42 bx = -32;
43 by = 48;
44 SCIPintervalSetBounds(&rhs, -21.333333336466665741681936197, -21.3333333341333322152877371991);
45 SCIPintervalSetBounds(&xbnds, -10.0382009139778674011722614523, 0.0);
46 SCIPintervalSetBounds(&ybnds, -7.93110393120117596055251851794, -0.537524704319278456843278490851);
47
48 // the quadratic equation
49 // ax*x^2 + ay*y^2 + axy*x*y + bx*x + by*y \in rhs
50 // does not have a solution for x \in xbnds and y \in ybnds
51 // however, relaxing rhs very slightly gives a solution
52 // SCIPintervalSolveBivariateQuadExpressionAllScalar() is expected to either
53 // - figure out that there is no solution, or
54 // - return compute some bounds because it decided to relax the bounds a bit
55 // However, it is not supposed to throw an assert (#1861).
56
57 SCIPintervalSolveBivariateQuadExpressionAllScalar(
58 infinity, /**< value for infinity in interval arithmetics */
59 &resultant, /**< buffer where to store result of operation */
60 ax, /**< square coefficient of x */
61 ay, /**< square coefficient of y */
62 axy, /**< bilinear coefficients */
63 bx, /**< linear coefficient of x */
64 by, /**< linear coefficient of y */
65 rhs, /**< right-hand-side of equation */
66 xbnds, /**< bounds on x */
67 ybnds /**< bounds on y */
68 );
69 }
70
71 /** some tests for SCIPintervalSolveUnivariateQuadExpression */
Test(intervalarith,solveuniquad)72 Test(intervalarith, solveuniquad)
73 {
74 SCIP_INTERVAL resultant;
75 SCIP_INTERVAL expect;
76 SCIP_INTERVAL sqrcoef;
77 SCIP_INTERVAL lincoef;
78 SCIP_INTERVAL rhs;
79 SCIP_INTERVAL xbnds;
80
81 SCIPintervalSetEntire(SCIP_DEFAULT_INFINITY, &xbnds);
82
83 // test with lincoef = 0
84 SCIPintervalSet(&lincoef, 0.0);
85 for( sqrcoef.inf = -2.0; sqrcoef.inf <= 2.0; ++sqrcoef.inf )
86 {
87 sqrcoef.sup = sqrcoef.inf;
88
89 for( rhs.inf = -2.0; rhs.inf <= 2.0; rhs.inf += 2.0 )
90 {
91 for( rhs.sup = rhs.inf + 1.0; rhs.sup <= 3.0; rhs.sup += 1.0 )
92 {
93 /* sqrcoef * x^2 = rhs -> x^2 = rhs/sqrcoef */
94 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
95 /* printf("%g x^2 in [%g,%g] -> [%g,%g]\n", sqrcoef.inf, rhs.inf, rhs.sup, resultant.inf, resultant.sup); */
96
97 if( sqrcoef.inf != 0.0 )
98 {
99 SCIPintervalDivScalar(SCIP_DEFAULT_INFINITY, &expect, rhs, sqrcoef.inf);
100 SCIPintervalSquareRoot(SCIP_DEFAULT_INFINITY, &expect, expect);
101 if( !SCIPintervalIsEmpty(SCIP_DEFAULT_INFINITY, expect) )
102 expect.inf = -expect.sup;
103 }
104 else if( rhs.inf <= 0.0 && rhs.sup >= 0.0 )
105 {
106 SCIPintervalSetEntire(SCIP_DEFAULT_INFINITY, &expect);
107 }
108 else
109 SCIPintervalSetEmpty(&expect);
110 /* printf("%g*x^2 = [%g,%g]; expect [%g,%g], got [%g,%g]\n", sqrcoef.inf, rhs.inf, rhs.sup, expect.inf, expect.sup, resultant.inf, resultant.sup); */
111
112 cr_assert(SCIPintervalIsEmpty(SCIP_DEFAULT_INFINITY, resultant) == SCIPintervalIsEmpty(SCIP_DEFAULT_INFINITY, expect));
113 cr_assert(SCIPintervalIsEntire(SCIP_DEFAULT_INFINITY, resultant) == SCIPintervalIsEntire(SCIP_DEFAULT_INFINITY, expect));
114 if( !SCIPintervalIsEmpty(SCIP_DEFAULT_INFINITY, expect) && !SCIPintervalIsEntire(SCIP_DEFAULT_INFINITY, expect) )
115 {
116 cr_assert_float_eq(resultant.inf, expect.inf, 1e-12, "unexpected x.inf %g for %g*x^2=[%g,%g], expected %g", resultant.inf, sqrcoef.inf, rhs.inf, rhs.sup, expect.inf);
117 cr_assert_float_eq(resultant.sup, expect.sup, 1e-12, "unexpected x.sup %g for %g*x^2=[%g,%g], expected %g", resultant.sup, sqrcoef.inf, rhs.inf, rhs.sup, expect.sup);
118 }
119 }
120 }
121 }
122
123 // test with lincoef != 0
124 for( lincoef.inf = -6.0; lincoef.inf <= 6.0; lincoef.inf += 3.0 )
125 {
126 lincoef.sup = lincoef.inf;
127
128 for( rhs.sup = -2.0; rhs.sup <= 2.0; rhs.sup += 2.0 )
129 {
130 rhs.inf = rhs.sup-1.0;
131 SCIPintervalSet(&sqrcoef, 0.5);
132 /* 0.5*x^2 + lincoef * x = rhs -> x = -lincoef +- sqrt(2*rhs + lincoef^2)
133 * thus: empty for 2*rhs+lincoef^2 < 0
134 */
135
136 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
137
138 if( 2.0*rhs.sup + lincoef.inf * lincoef.inf < 0.0 )
139 cr_assert(SCIPintervalIsEmpty(SCIP_DEFAULT_INFINITY, resultant));
140 else
141 {
142 expect.inf = -lincoef.inf - sqrt(2.0*rhs.sup + lincoef.inf*lincoef.inf);
143 expect.sup = -lincoef.inf + sqrt(2.0*rhs.sup + lincoef.inf*lincoef.inf);
144 cr_assert_float_eq(resultant.inf, expect.inf, 1e-12, "unexpected x.inf %g for 0.5*x^2%+g*x=[%g,%g], expected %g", resultant.inf, lincoef.inf, rhs.inf, rhs.sup, expect.inf);
145 cr_assert_float_eq(resultant.sup, expect.sup, 1e-12, "unexpected x.sup %g for 0.5*x^2%+g*x=[%g,%g], expected %g", resultant.sup, lincoef.inf, rhs.inf, rhs.sup, expect.sup);
146 }
147
148 SCIPintervalSet(&sqrcoef, -0.5);
149 /* -0.5*x^2 + lincoef * x = rhs -> x = lincoef +- sqrt(-2*rhs + lincoef^2)
150 * thus: empty for -2*rhs+lincoef^2 < 0
151 */
152
153 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
154
155 if( -2.0*rhs.inf + lincoef.inf * lincoef.inf < 0.0 )
156 cr_assert(SCIPintervalIsEmpty(SCIP_DEFAULT_INFINITY, resultant));
157 else
158 {
159 expect.inf = lincoef.inf - sqrt(-2.0*rhs.inf + lincoef.inf*lincoef.inf);
160 expect.sup = lincoef.inf + sqrt(-2.0*rhs.inf + lincoef.inf*lincoef.inf);
161 cr_assert_float_eq(resultant.inf, expect.inf, 1e-12, "unexpected x.inf %g for -0.5*x^2%+g*x=[%g,%g], expected %g", resultant.inf, lincoef.inf, rhs.inf, rhs.sup, expect.inf);
162 cr_assert_float_eq(resultant.sup, expect.sup, 1e-12, "unexpected x.sup %g for -0.5*x^2%+g*x=[%g,%g], expected %g", resultant.sup, lincoef.inf, rhs.inf, rhs.sup, expect.sup);
163 }
164 }
165 }
166
167 /* 0.5x^2-x has a minimum at x=1 (value -0.5) and maxima at the bounds of x
168 * further, {x : 0.5x^2-x >= 0} = [-infty,0] v [2,infty]
169 * further, {x : 0.5x^2-x in [0,1.5]} = [-1,0] v [2,3]
170 * further, {x : 0.5x^2-x >= 1.5} = [-infty,-1] v [3,infty]
171 */
172 SCIPintervalSet(&sqrcoef, 0.5);
173 SCIPintervalSet(&lincoef, -1.0);
174 SCIPintervalSetBounds(&rhs, 0.0, SCIP_DEFAULT_INFINITY);
175 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
176 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
177 cr_assert_eq(resultant.sup, +SCIP_DEFAULT_INFINITY);
178
179 SCIPintervalSetBounds(&rhs, 0.0, 1.5);
180 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
181 cr_assert_float_eq(resultant.inf, -1.0, 1e-12);
182 cr_assert_float_eq(resultant.sup, 3.0, 1e-12);
183
184 SCIPintervalSetBounds(&rhs, 1.5, SCIP_DEFAULT_INFINITY);
185 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
186 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
187 cr_assert_eq(resultant.sup, +SCIP_DEFAULT_INFINITY);
188
189 /* now, let's look only for solutions x >= 0:
190 * {x >= 0 : 0.5x^2-x >= 0} = [0,0] v [2,infty]
191 * {x >= 0 : 0.5x^2-x in [0,1.5]} = [0,0] v [2,3]
192 * {x >= 0 : 0.5x^2-x >= 1.5} = [3,infty]
193 */
194 SCIPintervalSetBounds(&xbnds, 0.0, SCIP_DEFAULT_INFINITY);
195
196 SCIPintervalSetBounds(&rhs, 0.0, SCIP_DEFAULT_INFINITY);
197 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
198 cr_assert_float_eq(resultant.inf, 0.0, 1e-12);
199 cr_assert_eq(resultant.sup, +SCIP_DEFAULT_INFINITY);
200
201 SCIPintervalSetBounds(&rhs, 0.0, 1.5);
202 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
203 cr_assert_float_eq(resultant.inf, 0.0, 1e-12);
204 cr_assert_float_eq(resultant.sup, 3.0, 1e-12);
205
206 SCIPintervalSetBounds(&rhs, 1.5, SCIP_DEFAULT_INFINITY);
207 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
208 cr_assert_float_eq(resultant.inf, 3.0, 1e-12);
209 cr_assert_eq(resultant.sup, +SCIP_DEFAULT_INFINITY);
210
211
212 /* now, let's look only for solutions x >= 1:
213 * {x >= 1 : 0.5x^2-x >= 0} = [2,infty]
214 * {x >= 1 : 0.5x^2-x in [0,1.5]} = [2,3]
215 * {x >= 1 : 0.5x^2-x >= 1.5} = [3,infty]
216 */
217 SCIPintervalSetBounds(&xbnds, 1.0, SCIP_DEFAULT_INFINITY);
218
219 SCIPintervalSetBounds(&rhs, 0.0, SCIP_DEFAULT_INFINITY);
220 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
221 cr_assert_float_eq(resultant.inf, 2.0, 1e-12);
222 cr_assert_eq(resultant.sup, +SCIP_DEFAULT_INFINITY);
223
224 SCIPintervalSetBounds(&rhs, 0.0, 1.5);
225 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
226 cr_assert_float_eq(resultant.inf, 2.0, 1e-12);
227 cr_assert_float_eq(resultant.sup, 3.0, 1e-12);
228
229 SCIPintervalSetBounds(&rhs, 1.5, SCIP_DEFAULT_INFINITY);
230 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
231 cr_assert_float_eq(resultant.inf, 3.0, 1e-12);
232 cr_assert_eq(resultant.sup, +SCIP_DEFAULT_INFINITY);
233
234
235 /* similarly, we can look only for solutions x <= -1:
236 * {x <= -1 : 0.5x^2-x >= 0} = [-infty,-1]
237 * {x <= -1 : 0.5x^2-x in [0,1.5]} = [-1,-1]
238 * {x <= -1 : 0.5x^2-x >= 1.5} = [-infty,-1]
239 */
240 SCIPintervalSetBounds(&xbnds, -SCIP_DEFAULT_INFINITY, -1.0);
241
242 SCIPintervalSetBounds(&rhs, 0.0, SCIP_DEFAULT_INFINITY);
243 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
244 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
245 cr_assert_float_eq(resultant.sup, -1.0, 1e-12);
246
247 SCIPintervalSetBounds(&rhs, 0.0, 1.5);
248 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
249 cr_assert_float_eq(resultant.inf, -1.0, 1e-12);
250 cr_assert_float_eq(resultant.sup, -1.0, 1e-12);
251
252 SCIPintervalSetBounds(&rhs, 1.5, SCIP_DEFAULT_INFINITY);
253 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
254 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
255 cr_assert_float_eq(resultant.sup, -1.0, 1e-12);
256
257
258 /* a linear equation with interval coefficients should be solved well:
259 * {x : [-1,1]*x >= 1.0} = [-infty,-1] v [1,infty]
260 * {x : [-1,1]*x <= 1.0} = [-infty,infty]
261 * {x : [-1,1]*x <=-1.0} = [-infty,-1] v [1,infty]
262 * {x : [-1,1]*x >=-1.0} = [-infty,infty]
263 *
264 * {x : [0,1]*x >= 1.0} = [1,infty]
265 * {x : [-1,0]*x >= 1.0} = [-infty,-1]
266 * {x : [0,1]*x >= 0.0} = [-infty,infty]
267 */
268 SCIPintervalSetEntire(SCIP_DEFAULT_INFINITY, &xbnds);
269 SCIPintervalSet(&sqrcoef, 0.0);
270 SCIPintervalSetBounds(&lincoef, -1.0, 1.0);
271 SCIPintervalSetBounds(&rhs, 1.0, SCIP_DEFAULT_INFINITY);
272 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
273 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
274 cr_assert_eq(resultant.sup, SCIP_DEFAULT_INFINITY);
275
276 SCIPintervalSetBounds(&rhs, -SCIP_DEFAULT_INFINITY, 1.0);
277 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
278 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
279 cr_assert_eq(resultant.sup, SCIP_DEFAULT_INFINITY);
280
281 SCIPintervalSetBounds(&rhs, -SCIP_DEFAULT_INFINITY, -1.0);
282 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
283 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
284 cr_assert_eq(resultant.sup, SCIP_DEFAULT_INFINITY);
285
286 SCIPintervalSetBounds(&rhs, -1.0, SCIP_DEFAULT_INFINITY);
287 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
288 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
289 cr_assert_eq(resultant.sup, SCIP_DEFAULT_INFINITY);
290
291 SCIPintervalSetBounds(&lincoef, 0.0, 1.0);
292 SCIPintervalSetBounds(&rhs, 1.0, SCIP_DEFAULT_INFINITY);
293 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
294 cr_assert_float_eq(resultant.inf, 1.0, 1e-12);
295 cr_assert_eq(resultant.sup, SCIP_DEFAULT_INFINITY);
296
297 SCIPintervalSetBounds(&lincoef, -1.0, 0.0);
298 SCIPintervalSetBounds(&rhs, 1.0, SCIP_DEFAULT_INFINITY);
299 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
300 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
301 cr_assert_float_eq(resultant.sup, -1.0, 1e-12);
302
303 SCIPintervalSetBounds(&lincoef, 0.0, 1.0);
304 SCIPintervalSetBounds(&rhs, 0.0, SCIP_DEFAULT_INFINITY);
305 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
306 cr_assert(SCIPintervalIsEntire(SCIP_DEFAULT_INFINITY, resultant));
307
308 /* providing bounds on x should result in better results:
309 * {x in [-infty,0]: [-1,1]*x >= 1.0} = [-infty,-1]
310 * {x in [0,+infty]: [-1,1]*x >= 1.0} = [1,infty]
311 * {x in [-infty,0]: [-1,1]*x <=-1.0} = [-infty,-1]
312 * {x in [0,+infty]: [-1,1]*x <=-1.0} = [1,infty]
313 */
314
315 SCIPintervalSetBounds(&lincoef, -1.0, 1.0);
316 SCIPintervalSetBounds(&rhs, 1.0, SCIP_DEFAULT_INFINITY);
317 SCIPintervalSetBounds(&xbnds, -SCIP_DEFAULT_INFINITY, 0.0);
318 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
319 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
320 cr_assert_float_eq(resultant.sup, -1.0, 1e-12);
321
322 SCIPintervalSetBounds(&xbnds, 0.0, SCIP_DEFAULT_INFINITY);
323 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
324 cr_assert_float_eq(resultant.inf, 1.0, 1e-12);
325 cr_assert_eq(resultant.sup, SCIP_DEFAULT_INFINITY);
326
327 SCIPintervalSetBounds(&rhs, -SCIP_DEFAULT_INFINITY, -1.0);
328 SCIPintervalSetBounds(&xbnds, -SCIP_DEFAULT_INFINITY, 0.0);
329 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
330 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
331 cr_assert_float_eq(resultant.sup, -1.0, 1e-12);
332
333 SCIPintervalSetBounds(&xbnds, 0.0, SCIP_DEFAULT_INFINITY);
334 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
335 cr_assert_float_eq(resultant.inf, 1.0, 1e-12);
336 cr_assert_eq(resultant.sup, SCIP_DEFAULT_INFINITY);
337
338 /* some more tests with lincoef=0 */
339 /* [0,1]*x^2 = 1 -> x^2 = [1,infty] -> x = [-infty,-1] v [1,infty]*/
340 SCIPintervalSetBounds(&sqrcoef, 0.0, 1.0);
341 SCIPintervalSet(&lincoef, 0.0);
342 SCIPintervalSetBounds(&rhs, 1.0, 1.0);
343 SCIPintervalSetEntire(SCIP_DEFAULT_INFINITY, &xbnds);
344 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
345 cr_assert(SCIPintervalIsEntire(SCIP_DEFAULT_INFINITY, resultant));
346
347 /* [0,1]*x^2 = 1, x >= 0 -> x >= 1 */
348 SCIPintervalSetBounds(&sqrcoef, 0.0, 1.0);
349 SCIPintervalSet(&lincoef, 0.0);
350 SCIPintervalSetBounds(&rhs, 1.0, 1.0);
351 SCIPintervalSetBounds(&xbnds, 0.0, SCIP_DEFAULT_INFINITY);
352 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
353 cr_assert_float_eq(resultant.inf, 1.0, 1e-12);
354 cr_assert_eq(resultant.sup, SCIP_DEFAULT_INFINITY);
355
356 /* [0,1]*x^2 = 1, x <= 0 -> x <= -1 */
357 SCIPintervalSetBounds(&sqrcoef, 0.0, 1.0);
358 SCIPintervalSet(&lincoef, 0.0);
359 SCIPintervalSetBounds(&rhs, 1.0, 1.0);
360 SCIPintervalSetBounds(&xbnds, -SCIP_DEFAULT_INFINITY, -1.0);
361 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, lincoef, rhs, xbnds);
362 cr_assert_eq(resultant.inf, -SCIP_DEFAULT_INFINITY);
363 cr_assert_float_eq(resultant.sup, -1.0, 1e-12);
364 }
365
366 /** some tests for x*y in rhs */
Test(intervalarith,xy)367 Test(intervalarith, xy)
368 {
369 SCIP_INTERVAL xbnds;
370 SCIP_INTERVAL ybnds;
371 SCIP_INTERVAL rhs;
372 SCIP_INTERVAL resultant;
373 SCIP_INTERVAL sqrcoef;
374
375 /* x*y == 1.0 for x,y in [-1,1] -> x = [-1,-1] v [1,1] */
376 SCIPintervalSetBounds(&xbnds, -1.0, 1.0);
377 SCIPintervalSetBounds(&ybnds, -1.0, 1.0);
378 SCIPintervalSetBounds(&rhs, 1.0, 1.0);
379 SCIPintervalSolveBivariateQuadExpressionAllScalar(SCIP_DEFAULT_INFINITY, &resultant, 0.0, 0.0, 1.0, 0.0, 0.0, rhs, xbnds, ybnds);
380 cr_assert_eq(resultant.inf, -1.0);
381 cr_assert_eq(resultant.sup, 1.0);
382
383 /* with x in [-1,0], this should then give x = [-1,-1] */
384 SCIPintervalSetBounds(&xbnds, -1.0, 0.0);
385 SCIPintervalSolveBivariateQuadExpressionAllScalar(SCIP_DEFAULT_INFINITY, &resultant, 0.0, 0.0, 1.0, 0.0, 0.0, rhs, xbnds, ybnds);
386 cr_assert_eq(resultant.inf, -1.0);
387 cr_assert_float_eq(resultant.sup, -1.0, 1e-12);
388
389 /* and with x in [0,1], this should then give x = [1,1] */
390 SCIPintervalSetBounds(&xbnds, 0.0, 1.0);
391 SCIPintervalSolveBivariateQuadExpressionAllScalar(SCIP_DEFAULT_INFINITY, &resultant, 0.0, 0.0, 1.0, 0.0, 0.0, rhs, xbnds, ybnds);
392 cr_assert_float_eq(resultant.inf, 1.0, 1e-12);
393 cr_assert_eq(resultant.sup, 1.0);
394
395 /* x*y >= 1.0 for x,y in [-1,1] -> x = [-1,-1] */
396 SCIPintervalSetBounds(&xbnds, -1.0, 1.0);
397 SCIPintervalSetBounds(&ybnds, -1.0, 1.0);
398 SCIPintervalSetBounds(&rhs, 1.0, SCIP_DEFAULT_INFINITY);
399 SCIPintervalSolveBivariateQuadExpressionAllScalar(SCIP_DEFAULT_INFINITY, &resultant, 0.0, 0.0, 1.0, 0.0, 0.0, rhs, xbnds, ybnds);
400 cr_assert_eq(resultant.inf, -1.0);
401 cr_assert_float_eq(resultant.sup, 1.0, 1e-12);
402
403 /* x*y >= 1.0 for x in [-1,1], y in [-1,0] -> x = [-1,-1] */
404 SCIPintervalSetBounds(&ybnds, -1.0, -0.0);
405 SCIPintervalSolveBivariateQuadExpressionAllScalar(SCIP_DEFAULT_INFINITY, &resultant, 0.0, 0.0, 1.0, 0.0, 0.0, rhs, xbnds, ybnds);
406 cr_assert_eq(resultant.inf, -1.0);
407 /* currently still gives 1 as upper bound, so cr_assert_float_eq(resultant.sup, -1.0, 1e-12); fails
408 * however, SCIPintervalSolveUnivariateQuadExpression handles this better:
409 */
410 SCIPintervalSet(&sqrcoef, 0.0);
411 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, ybnds, rhs, xbnds);
412 cr_assert_eq(resultant.inf, -1.0);
413 cr_assert_float_eq(resultant.sup, -1.0, 1e-12);
414
415 /* similar for x*y >= 1.0 for x in [-1,1], y in [0,1] -> x = [1,1] */
416 SCIPintervalSetBounds(&ybnds, 0, 1.0);
417 SCIPintervalSolveBivariateQuadExpressionAllScalar(SCIP_DEFAULT_INFINITY, &resultant, 0.0, 0.0, 1.0, 0.0, 0.0, rhs, xbnds, ybnds);
418 cr_assert_eq(resultant.sup, 1.0);
419 /* currently still gives -1 as lower bound, so cr_assert_float_eq(resultant.inf, 1.0, 1e-12); fails
420 * however, SCIPintervalSolveUnivariateQuadExpression handles this better:
421 */
422 SCIPintervalSet(&sqrcoef, 0.0);
423 SCIPintervalSolveUnivariateQuadExpression(SCIP_DEFAULT_INFINITY, &resultant, sqrcoef, ybnds, rhs, xbnds);
424 cr_assert_float_eq(resultant.inf, 1.0, 1e-12);
425 cr_assert_eq(resultant.sup, 1.0);
426 }
427
Test(intervalarith,issue2250)428 Test(intervalarith, issue2250)
429 {
430 SCIP_Real infinity;
431 SCIP_INTERVAL resultant;
432 SCIP_Real ax;
433 SCIP_Real ay;
434 SCIP_Real axy;
435 SCIP_Real bx;
436 SCIP_Real by;
437 SCIP_INTERVAL rhs;
438 SCIP_INTERVAL xbnds;
439 SCIP_INTERVAL ybnds;
440
441 infinity = 1e43;
442 ax = 1;
443 ay = 0;
444 axy = 1;
445 bx = -2;
446 by = -6;
447 SCIPintervalSetBounds(&rhs, -infinity, 0.0);
448 SCIPintervalSetBounds(&xbnds, 0.0, infinity);
449 SCIPintervalSetBounds(&ybnds, 0.0, infinity);
450
451 /* x=y=4 is feasible for this equation, so x=4 should be part of the solution of this equation */
452 cr_assert(ax*4*4 + bx*4 + ay*4*4 + by*4 + axy*4*4 <= 1e-12);
453
454 SCIPintervalSolveBivariateQuadExpressionAllScalar(
455 infinity, /**< value for infinity in interval arithmetics */
456 &resultant, /**< buffer where to store result of operation */
457 ax, /**< square coefficient of x */
458 ay, /**< square coefficient of y */
459 axy, /**< bilinear coefficients */
460 bx, /**< linear coefficient of x */
461 by, /**< linear coefficient of y */
462 rhs, /**< right-hand-side of equation */
463 xbnds, /**< bounds on x */
464 ybnds /**< bounds on y */
465 );
466
467 cr_assert(resultant.inf <= 4.0);
468 cr_assert(resultant.sup >= 4.0);
469 }
470
471 /* The fail in #2650 was caused by GCC reorganizing operations in SCIPintervalReciprocal so that divisions
472 * were not evaluated with the correct rounding mode.
473 * Unfortunately, I was not able to reproduce this with a single test like this
474 * (build SCIP with OPT=dbg SHARED=true USRCFLAGS="-O3 -DNDEBUG -fomit-frame-pointer").
475 */
Test(intervalarith,issue2650)476 Test(intervalarith, issue2650)
477 {
478 SCIP_Real infinity = 1.0e300;
479 SCIP_INTERVAL resultant;
480 SCIP_INTERVAL operand;
481 SCIP_INTERVAL base;
482 SCIP_INTERVAL image;
483
484 SCIPintervalSetBounds(&base, 0.0, infinity);
485 SCIPintervalSetBounds(&image, 9.0, 81.0);
486 SCIPintervalPowerScalarInverse(infinity, &resultant, base, 0.2, image);
487 printf("x^0.2 = [%.15g,%.15g] -> x = [%.15g,%.15g]\n", image.inf, image.sup, resultant.inf, resultant.sup);
488
489 cr_assert(resultant.inf <= 3486784401.0);
490 cr_assert(resultant.sup >= 3486784401.0);
491
492 /* the code above failed because the 1/0.2 wasn't computed correctly: */
493 SCIPintervalSetBounds(&operand, 0.2, 0.2);
494 SCIPintervalReciprocal(infinity, &resultant, operand);
495 printf("1/[0.2,0.2] = [%.15g,%.15g]\n", resultant.inf, resultant.sup);
496
497 cr_assert(resultant.inf <= 5.0);
498 cr_assert(resultant.sup >= 5.0);
499 }
500