1 /**
2 * Copyright 2015, SRI International.
3 *
4 * This file is part of LibPoly.
5 *
6 * LibPoly is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * LibPoly is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with LibPoly. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 #include <interval.h>
21
22 #include "interval/arithmetic.h"
23
24 #include "number/rational.h"
25 #include "number/dyadic_rational.h"
26 #include "number/value.h"
27
rational_interval_endpoint_lt(const lp_rational_t * a,int a_open,const lp_rational_t * b,int b_open)28 int rational_interval_endpoint_lt(const lp_rational_t* a, int a_open, const lp_rational_t* b, int b_open) {
29 int cmp = rational_cmp(a, b);
30 if (cmp == 0) {
31 return (!a_open && b_open);
32 } else {
33 return cmp < 0;
34 }
35 }
36
dyadic_interval_endpoint_lt(const lp_dyadic_rational_t * a,int a_open,const lp_dyadic_rational_t * b,int b_open)37 int dyadic_interval_endpoint_lt(const lp_dyadic_rational_t* a, int a_open, const lp_dyadic_rational_t*b, int b_open) {
38 int cmp = dyadic_rational_cmp(a, b);
39 if (cmp == 0) {
40 return (!a_open && b_open);
41 } else {
42 return cmp < 0;
43 }
44 }
45
lp_interval_endpoint_lt(const lp_value_t * a,int a_open,const lp_value_t * b,int b_open)46 int lp_interval_endpoint_lt(const lp_value_t* a, int a_open, const lp_value_t* b, int b_open) {
47 int cmp = lp_value_cmp(a, b);
48 if (cmp == 0) {
49 return (!a_open && b_open);
50 } else {
51 return cmp < 0;
52 }
53 }
54
rational_interval_add(lp_rational_interval_t * S,const lp_rational_interval_t * I1,const lp_rational_interval_t * I2)55 void rational_interval_add(lp_rational_interval_t* S, const lp_rational_interval_t* I1, const lp_rational_interval_t* I2) {
56
57 if (I1->is_point && I2->is_point) {
58 if (!S->is_point) {
59 rational_destruct(&S->b);
60 }
61 rational_add(&S->a, &I1->a, &I2->a);
62 S->b_open = S->a_open = 0;
63 S->is_point = 1;
64 return;
65 }
66
67 if (I2->is_point) {
68 // Reuse symmetry
69 rational_interval_add(S, I2, I1);
70 return;
71 }
72
73 lp_rational_interval_t result;
74
75 if (I1->is_point) {
76 // Just shift by I1->a (I2 is not a point)
77 lp_rational_interval_construct_copy(&result, I2);
78 rational_add(&result.a, &result.a, &I1->a);
79 rational_add(&result.b, &result.b, &I1->a);
80 } else {
81 // [a, b] + [c, d] = [a + c, b + d]
82 rational_construct(&result.a);
83 rational_construct(&result.b);
84 rational_add(&result.a, &I1->a, &I2->a);
85 rational_add(&result.b, &I1->b, &I2->b);
86 result.a_open = I1->a_open || I2->a_open;
87 result.b_open = I1->b_open || I2->b_open;
88 result.is_point = 0;
89 }
90
91 lp_rational_interval_swap(&result, S);
92 lp_rational_interval_destruct(&result);
93 }
94
dyadic_interval_add(lp_dyadic_interval_t * S,const lp_dyadic_interval_t * I1,const lp_dyadic_interval_t * I2)95 void dyadic_interval_add(lp_dyadic_interval_t* S, const lp_dyadic_interval_t* I1, const lp_dyadic_interval_t* I2) {
96
97 if (I1->is_point && I2->is_point) {
98 if (!S->is_point) {
99 dyadic_rational_destruct(&S->b);
100 }
101 dyadic_rational_add(&S->a, &I1->a, &I2->a);
102 S->b_open = S->a_open = 0;
103 S->is_point = 1;
104 return;
105 }
106
107 if (I2->is_point) {
108 // Reuse symmetry
109 dyadic_interval_add(S, I2, I1);
110 return;
111 }
112
113 lp_dyadic_interval_t result;
114
115 if (I1->is_point) {
116 // Just shift by I1->a
117 lp_dyadic_interval_construct_copy(&result, I2);
118 dyadic_rational_add(&result.a, &result.a, &I1->a);
119 dyadic_rational_add(&result.b, &result.b, &I1->a);
120 } else {
121 // Both non-points
122 // [a, b] + [c, d] = [a + c, b + d]
123 dyadic_rational_construct(&result.a);
124 dyadic_rational_construct(&result.b);
125 dyadic_rational_add(&result.a, &I1->a, &I2->a);
126 dyadic_rational_add(&result.b, &I1->b, &I2->b);
127 result.a_open = I1->a_open || I2->a_open;
128 result.b_open = I1->b_open || I2->b_open;
129 result.is_point = 0;
130 }
131
132 lp_dyadic_interval_swap(&result, S);
133 lp_dyadic_interval_destruct(&result);
134 }
135
136 static
lp_value_add_approx(const lp_value_t * v1,const lp_value_t * v2,lp_value_t * lb,lp_value_t * ub)137 int lp_value_add_approx(const lp_value_t* v1, const lp_value_t* v2, lp_value_t* lb, lp_value_t* ub) {
138
139 assert(v1 != lb && v1 != ub);
140 assert(v2 != lb && v2 != ub);
141 assert(lb != ub);
142
143 lp_integer_t add_int;
144 lp_rational_t add_rat;
145 lp_dyadic_rational_t add_dy;
146 lp_dyadic_interval_t add_dy_interval;
147
148 int is_point = 1;
149
150 assert(v1->type != LP_VALUE_NONE);
151 assert(v2->type != LP_VALUE_NONE);
152
153 if (v1->type == v2->type) {
154 switch (v1->type) {
155 case LP_VALUE_MINUS_INFINITY:
156 if (lb) {
157 lp_value_assign_raw(lb, LP_VALUE_MINUS_INFINITY, 0);
158 }
159 if (ub) {
160 lp_value_assign_raw(ub, LP_VALUE_MINUS_INFINITY, 0);
161 }
162 break;
163 case LP_VALUE_PLUS_INFINITY:
164 if (lb) {
165 lp_value_assign_raw(lb, LP_VALUE_PLUS_INFINITY, 0);
166 }
167 if (ub) {
168 lp_value_assign_raw(ub, LP_VALUE_PLUS_INFINITY, 0);
169 }
170 break;
171 case LP_VALUE_INTEGER:
172 lp_integer_construct(&add_int);
173 lp_integer_add(lp_Z, &add_int, &v1->value.z, &v2->value.z);
174 if (lb) {
175 lp_value_assign_raw(lb, LP_VALUE_INTEGER, &add_int);
176 }
177 if (ub) {
178 lp_value_assign_raw(ub, LP_VALUE_INTEGER, &add_int);
179 }
180 lp_integer_destruct(&add_int);
181 break;
182 case LP_VALUE_RATIONAL:
183 lp_rational_construct(&add_rat);
184 lp_rational_add(&add_rat, &v1->value.q, &v2->value.q);
185 if (lb) {
186 lp_value_assign_raw(lb, LP_VALUE_RATIONAL, &add_rat);
187 }
188 if (ub) {
189 lp_value_assign_raw(ub, LP_VALUE_RATIONAL, &add_rat);
190 }
191 lp_rational_destruct(&add_rat);
192 break;
193 case LP_VALUE_DYADIC_RATIONAL:
194 lp_dyadic_rational_construct(&add_dy);
195 lp_dyadic_rational_add(&add_dy, &v1->value.dy_q, &v2->value.dy_q);
196 if (lb) {
197 lp_value_assign_raw(lb, LP_VALUE_DYADIC_RATIONAL, &add_dy);
198 }
199 if (ub) {
200 lp_value_assign_raw(ub, LP_VALUE_DYADIC_RATIONAL, &add_dy);
201 }
202 lp_dyadic_rational_destruct(&add_dy);
203 break;
204 case LP_VALUE_ALGEBRAIC:
205 if (v1->value.a.I.is_point && v2->value.a.I.is_point) {
206 lp_dyadic_rational_construct(&add_dy);
207 lp_dyadic_rational_add(&add_dy, &v1->value.a.I.a, &v2->value.a.I.a);
208 if (lb) {
209 lp_value_assign_raw(lb, LP_VALUE_DYADIC_RATIONAL, &add_dy);
210 }
211 if (ub) {
212 lp_value_assign_raw(ub, LP_VALUE_DYADIC_RATIONAL, &add_dy);
213 }
214 lp_dyadic_rational_destruct(&add_dy);
215 } else {
216 lp_dyadic_interval_construct_zero(&add_dy_interval);
217 dyadic_interval_add(&add_dy_interval, &v1->value.a.I, &v2->value.a.I);
218 if (lb) {
219 lp_value_assign_raw(lb, LP_VALUE_DYADIC_RATIONAL, &add_dy_interval.a);
220 }
221 if (ub) {
222 lp_value_assign_raw(ub, LP_VALUE_DYADIC_RATIONAL, &add_dy_interval.b);
223 }
224 is_point = 0;
225 }
226 break;
227 case LP_VALUE_NONE:
228 assert(0);
229 break;
230 }
231 } else {
232
233 lp_value_t v1_tmp, v2_tmp;
234 const lp_value_t* v1_to_use;
235 const lp_value_t* v2_to_use;
236
237 // Cast to the same type
238 int cast = lp_value_to_same_type(v1, v2, &v1_tmp, &v2_tmp, &v1_to_use, &v2_to_use);
239 if (cast) {
240 // Cast successful, just add as usual
241 lp_value_add_approx(v1_to_use, v2_to_use, lb, ub);
242 // If a fesh value was used, delete it
243 if (v1_to_use != v1) {
244 lp_value_destruct(&v1_tmp);
245 }
246 if (v2_to_use != v2) {
247 lp_value_destruct(&v2_tmp);
248 }
249 } else {
250 // Cast failed, check the cases
251 int undefined = 0;
252 lp_value_type_t result = LP_VALUE_NONE;
253
254 // Check infinity cases
255 if (v1->type == LP_VALUE_MINUS_INFINITY) {
256 if (v2->type != LP_VALUE_PLUS_INFINITY) {
257 result = LP_VALUE_MINUS_INFINITY;
258 } else {
259 undefined = 1;
260 }
261 }
262 if (v1->type == LP_VALUE_PLUS_INFINITY) {
263 if (v2->type != LP_VALUE_MINUS_INFINITY) {
264 result = LP_VALUE_PLUS_INFINITY;
265 } else {
266 undefined = 1;
267 }
268 }
269 if (v2->type == LP_VALUE_MINUS_INFINITY) {
270 if (v1->type != LP_VALUE_PLUS_INFINITY) {
271 result = LP_VALUE_MINUS_INFINITY;
272 } else {
273 undefined = 1;
274 }
275 }
276 if (v2->type == LP_VALUE_PLUS_INFINITY) {
277 if (v1->type != LP_VALUE_MINUS_INFINITY) {
278 result = LP_VALUE_PLUS_INFINITY;
279 } else {
280 undefined = 1;
281 }
282 }
283
284 // If infinity handled, we're done
285 if (result != LP_VALUE_NONE) {
286 if (lb) {
287 lp_value_assign_raw(lb, result, 0);
288 }
289 if (ub) {
290 lp_value_assign_raw(ub, result, 0);
291 }
292 return 1;
293 }
294 // If undefined, return NONE
295 if (undefined) {
296 if (lb) {
297 lp_value_assign_raw(lb, LP_VALUE_NONE, 0);
298 }
299 if (ub) {
300 lp_value_assign_raw(ub, LP_VALUE_NONE, 0);
301 }
302 return 0;
303 }
304
305 // OK, now deal with algebraic numbers, assume v1 is algebraic
306 if (v2->type == LP_VALUE_ALGEBRAIC) {
307 const lp_value_t* tmp = v1; v1 = v2; v2 = tmp;
308 }
309
310 assert(v1->type == LP_VALUE_ALGEBRAIC);
311 assert(v2->type != LP_VALUE_MINUS_INFINITY);
312 assert(v2->type != LP_VALUE_PLUS_INFINITY);
313 assert(v2->type != LP_VALUE_ALGEBRAIC);
314
315 // v1 is of the form (dy1, dy2), v2 is a proper non-algebraic value
316 // we add them into (dy1 + v2, dy2 + v2)
317
318 lp_value_t v1_lb, v1_ub;
319 const lp_dyadic_rational_t* a = &v1->value.a.I.a;
320 const lp_dyadic_rational_t* b = v1->value.a.I.is_point ? a : &v1->value.a.I.b;
321 lp_value_construct(&v1_lb, LP_VALUE_DYADIC_RATIONAL, a);
322 lp_value_construct(&v1_ub, LP_VALUE_DYADIC_RATIONAL, b);
323 lp_value_add_approx(&v1_lb, v2, lb, 0);
324 lp_value_add_approx(&v1_ub, v2, ub, 0);
325 is_point = 0;
326 lp_value_destruct(&v1_lb);
327 lp_value_destruct(&v1_ub);
328 }
329 }
330
331 return is_point;
332 }
333
lp_interval_add(lp_interval_t * add,const lp_interval_t * I1,const lp_interval_t * I2)334 void lp_interval_add(lp_interval_t* add, const lp_interval_t* I1, const lp_interval_t* I2) {
335
336 lp_interval_t result;
337 lp_interval_construct_full(&result);
338
339 if (I1->is_point && I2->is_point) {
340 result.is_point = lp_value_add_approx(&I1->a, &I2->a, &result.a, &result.b);
341 if (result.is_point) {
342 lp_value_destruct(&result.b);
343 }
344 result.b_open = result.a_open = !result.is_point;
345 } else {
346 // [a, b] + [c, d] = [a + c, b + d]
347 const lp_value_t* I1_a = &I1->a;
348 const lp_value_t* I1_b = I1->is_point ? I1_a : &I1->b;
349 const lp_value_t* I2_a = &I2->a;
350 const lp_value_t* I2_b = I2->is_point ? I2_a : &I2->b;
351 int a_point = lp_value_add_approx(I1_a, I2_a, &result.a, 0);
352 int b_point = lp_value_add_approx(I1_b, I2_b, 0, &result.b);
353 result.a_open = I1->a_open || I2->a_open || !a_point;
354 result.b_open = I1->b_open || I2->b_open || !b_point;
355 result.is_point = 0;
356 }
357
358 lp_interval_swap(add, &result);
359 lp_interval_destruct(&result);
360 }
361
rational_interval_neg(lp_rational_interval_t * N,const lp_rational_interval_t * I)362 void rational_interval_neg(lp_rational_interval_t* N, const lp_rational_interval_t* I) {
363 if (I->is_point) {
364 if (!N->is_point) {
365 rational_destruct(&N->b);
366 }
367 rational_neg(&N->a, &I->a);
368 N->b_open = N->a_open = 0;
369 N->is_point = 1;
370 return;
371 }
372
373 if (N->is_point) {
374 rational_construct(&N->b);
375 N->is_point = 0;
376 }
377
378 // -[a, b] = [-b, -a]
379 // (doing swap in case I and N are the same)
380 rational_neg(&N->a, &I->a);
381 rational_neg(&N->b, &I->b);
382 N->a_open = I->a_open;
383 N->b_open = I->b_open;
384
385 rational_swap(&N->a, &N->b);
386 size_t tmp = N->a_open;
387 N->a_open = N->b_open;
388 N->b_open = tmp;
389 }
390
dyadic_interval_neg(lp_dyadic_interval_t * N,const lp_dyadic_interval_t * I)391 void dyadic_interval_neg(lp_dyadic_interval_t* N, const lp_dyadic_interval_t* I) {
392 if (I->is_point) {
393 if (!N->is_point) {
394 dyadic_rational_destruct(&N->b);
395 }
396 dyadic_rational_neg(&N->a, &I->a);
397 N->b_open = N->a_open = 0;
398 N->is_point = 1;
399 return;
400 }
401
402 if (N->is_point) {
403 dyadic_rational_construct(&N->b);
404 N->is_point = 0;
405 }
406
407 // -[a, b] = [-b, -a]
408 // (doing swap in case I and N are the same)
409 dyadic_rational_neg(&N->a, &I->a);
410 dyadic_rational_neg(&N->b, &I->b);
411 N->a_open = I->a_open;
412 N->b_open = I->b_open;
413
414 dyadic_rational_swap(&N->a, &N->b);
415 size_t tmp = N->a_open;
416 N->a_open = N->b_open;
417 N->b_open = tmp;
418 }
419
rational_interval_sub(lp_rational_interval_t * S,const lp_rational_interval_t * I1,const lp_rational_interval_t * I2)420 void rational_interval_sub(lp_rational_interval_t* S, const lp_rational_interval_t* I1, const lp_rational_interval_t* I2) {
421 lp_rational_interval_t neg;
422 lp_rational_interval_construct_copy(&neg, I2);
423 rational_interval_neg(&neg, &neg);
424 rational_interval_add(S, I1, &neg);
425 lp_rational_interval_destruct(&neg);
426 }
427
dyadic_interval_sub(lp_dyadic_interval_t * S,const lp_dyadic_interval_t * I1,const lp_dyadic_interval_t * I2)428 void dyadic_interval_sub(lp_dyadic_interval_t* S, const lp_dyadic_interval_t* I1, const lp_dyadic_interval_t* I2) {
429 lp_dyadic_interval_t neg;
430 lp_dyadic_interval_construct_copy(&neg, I2);
431 dyadic_interval_neg(&neg, &neg);
432 dyadic_interval_add(S, I1, &neg);
433 lp_dyadic_interval_destruct(&neg);
434 }
435
rational_interval_mul(lp_rational_interval_t * P,const lp_rational_interval_t * I1,const lp_rational_interval_t * I2)436 void rational_interval_mul(lp_rational_interval_t* P, const lp_rational_interval_t* I1, const lp_rational_interval_t* I2) {
437 if (I1->is_point) {
438 if (I2->is_point) {
439 // Just multiply the points
440 rational_mul(&P->a, &I1->a, &I2->a);
441 if (!P->is_point) {
442 rational_destruct(&P->b);
443 P->is_point = 1;
444 }
445 P->a_open = P->b_open = 0;
446 } else {
447 // I1 is point
448 // I2 is not point
449
450 // Depending on the sign of a, we might have to flip
451 int a_sgn = rational_sgn(&I1->a);
452 if (a_sgn == 0) {
453 // It's just 0
454 if (!P->is_point) {
455 rational_destruct(&P->b);
456 P->is_point = 1;
457 }
458 P->a_open = P->b_open = 0;
459 rational_assign_int(&P->a, 0, 1);
460 } else if (a_sgn > 0) {
461 lp_rational_interval_t result;
462 rational_construct(&result.a);
463 rational_construct(&result.b);
464 result.is_point = 0;
465 result.a_open = I2->a_open;
466 result.b_open = I2->b_open;
467 rational_mul(&result.a, &I1->a, &I2->a);
468 rational_mul(&result.b, &I1->a, &I2->b);
469 lp_rational_interval_swap(&result, P);
470 lp_rational_interval_destruct(&result);
471 } else {
472 lp_rational_interval_t result;
473 rational_construct(&result.a);
474 rational_construct(&result.b);
475 result.is_point = 0;
476 result.a_open = I2->a_open;
477 result.b_open = I2->b_open;
478 rational_mul(&result.b, &I1->a, &I2->a);
479 rational_mul(&result.a, &I1->a, &I2->b);
480 lp_rational_interval_swap(&result, P);
481 lp_rational_interval_destruct(&result);
482 }
483 }
484 } else if (I2->is_point) {
485 rational_interval_mul(P, I2, I1);
486 } else {
487 if (P->is_point) {
488 rational_construct(&P->b);
489 P->is_point = 0;
490 }
491
492 //
493 // I1 x I2 = { x*y | x in I1, y in I2 }
494 // = { x*y | I1.a < x < I1.b, I2.a < y < I2.b }
495 //
496
497 lp_rational_interval_t result;
498 lp_rational_interval_construct_zero(&result);
499
500 lp_rational_t tmp;
501 rational_construct(&tmp);
502
503 // I1.a x I2.a
504 rational_mul(&result.a, &I1->a, &I2->a);
505 rational_construct_copy(&result.b, &result.a);
506 result.a_open = result.b_open = I1->a_open || I2->a_open;
507 result.is_point = 0;
508
509 // I1.a x I2.b
510 int tmp_open = I1->a_open || I2->b_open;
511 rational_mul(&tmp, &I1->a, &I2->b);
512 if (rational_interval_endpoint_lt(&tmp, tmp_open, &result.a, result.a_open)) {
513 rational_swap(&tmp, &result.a);
514 result.a_open = tmp_open;
515 } else if (rational_interval_endpoint_lt(&result.b, result.b_open, &tmp, tmp_open)) {
516 rational_swap(&tmp, &result.b);
517 result.b_open = tmp_open;
518 }
519
520 // I1.b x I2.a
521 tmp_open = I1->b_open || I2->a_open;
522 rational_mul(&tmp, &I1->b, &I2->a);
523 if (rational_interval_endpoint_lt(&tmp, tmp_open, &result.a, result.a_open)) {
524 rational_swap(&tmp, &result.a);
525 result.a_open = tmp_open;
526 } else if (rational_interval_endpoint_lt(&result.b, result.b_open, &tmp, tmp_open)) {
527 rational_swap(&tmp, &result.b);
528 result.b_open = tmp_open;
529 }
530
531 // I1.b x I2.b
532 tmp_open = I1->b_open || I2->b_open;
533 rational_mul(&tmp, &I1->b, &I2->b);
534 if (rational_interval_endpoint_lt(&tmp, tmp_open, &result.a, result.a_open)) {
535 rational_swap(&tmp, &result.a);
536 result.a_open = tmp_open;
537 } else if (rational_interval_endpoint_lt(&result.b, result.b_open, &tmp, tmp_open)) {
538 rational_swap(&tmp, &result.b);
539 result.b_open = tmp_open;
540 }
541
542 lp_rational_interval_swap(&result, P);
543 lp_rational_interval_destruct(&result);
544 rational_destruct(&tmp);
545 }
546 }
547
dyadic_interval_mul(lp_dyadic_interval_t * P,const lp_dyadic_interval_t * I1,const lp_dyadic_interval_t * I2)548 void dyadic_interval_mul(lp_dyadic_interval_t* P, const lp_dyadic_interval_t* I1, const lp_dyadic_interval_t* I2) {
549 if (I1->is_point) {
550 if (I2->is_point) {
551 // Just multiply the points
552 dyadic_rational_mul(&P->a, &I1->a, &I2->a);
553 if (!P->is_point) {
554 dyadic_rational_destruct(&P->b);
555 P->is_point = 1;
556 }
557 P->a_open = P->b_open = 0;
558 } else {
559 // Depending on the sign of a, we might have to flip
560 int a_sgn = dyadic_rational_sgn(&I1->a);
561 if (a_sgn == 0) {
562 // It's just 0
563 if (!P->is_point) {
564 dyadic_rational_destruct(&P->b);
565 P->is_point = 1;
566 }
567 P->a_open = P->b_open = 0;
568 dyadic_rational_assign_int(&P->a, 0, 1);
569 } else if (a_sgn > 0) {
570 // Regular multiplication
571 lp_dyadic_interval_t result;
572 dyadic_rational_construct(&result.a);
573 dyadic_rational_construct(&result.b);
574 result.is_point = 0;
575 result.a_open = I2->a_open;
576 result.b_open = I2->b_open;
577 dyadic_rational_mul(&result.a, &I1->a, &I2->a);
578 dyadic_rational_mul(&result.b, &I1->a, &I2->b);
579 lp_dyadic_interval_swap(&result, P);
580 lp_dyadic_interval_destruct(&result);
581 } else {
582 lp_dyadic_interval_t result;
583 dyadic_rational_construct(&result.a);
584 dyadic_rational_construct(&result.b);
585 result.is_point = 0;
586 result.a_open = I2->b_open;
587 result.b_open = I2->a_open;
588 dyadic_rational_mul(&result.a, &I1->a, &I2->b);
589 dyadic_rational_mul(&result.b, &I1->a, &I2->a);
590 lp_dyadic_interval_swap(&result, P);
591 lp_dyadic_interval_destruct(&result);
592 }
593 }
594 } else if (I2->is_point) {
595 dyadic_interval_mul(P, I2, I1);
596 } else {
597 if (P->is_point) {
598 dyadic_rational_construct(&P->b);
599 P->is_point = 0;
600 }
601
602 //
603 // I1 x I2 = { x*y | x in I1, y in I2 }
604 // = { x*y | I1.a < x < I1.b, I2.a < y < I2.b }
605 //
606
607 lp_dyadic_interval_t result;
608 lp_dyadic_interval_construct_zero(&result);
609
610 lp_dyadic_rational_t tmp;
611 dyadic_rational_construct(&tmp);
612
613 // I1.a x I2.a
614 dyadic_rational_mul(&result.a, &I1->a, &I2->a);
615 dyadic_rational_construct_copy(&result.b, &result.a); // was not constructed (is 0)
616 result.a_open = result.b_open = I1->a_open || I2->a_open;
617 result.is_point = 0;
618
619 // I1.a x I2.b
620 int tmp_open = I1->a_open || I2->b_open;
621 dyadic_rational_mul(&tmp, &I1->a, &I2->b);
622 if (dyadic_interval_endpoint_lt(&tmp, tmp_open, &result.a, result.a_open)) {
623 dyadic_rational_swap(&tmp, &result.a);
624 result.a_open = tmp_open;
625 } else if (dyadic_interval_endpoint_lt(&result.b, result.b_open, &tmp, tmp_open)) {
626 dyadic_rational_swap(&tmp, &result.b);
627 result.b_open = tmp_open;
628 }
629
630 // I1.b x I2.a
631 tmp_open = I1->b_open || I2->a_open;
632 dyadic_rational_mul(&tmp, &I1->b, &I2->a);
633 if (dyadic_interval_endpoint_lt(&tmp, tmp_open, &result.a, result.a_open)) {
634 dyadic_rational_swap(&tmp, &result.a);
635 result.a_open = tmp_open;
636 } else if (dyadic_interval_endpoint_lt(&result.b, result.b_open, &tmp, tmp_open)) {
637 dyadic_rational_swap(&tmp, &result.b);
638 result.b_open = tmp_open;
639 }
640
641 // I1.b x I2.b
642 tmp_open = I1->b_open || I2->b_open;
643 dyadic_rational_mul(&tmp, &I1->b, &I2->b);
644 if (dyadic_interval_endpoint_lt(&tmp, tmp_open, &result.a, result.a_open)) {
645 dyadic_rational_swap(&tmp, &result.a);
646 result.a_open = tmp_open;
647 } else if (dyadic_interval_endpoint_lt(&result.b, result.b_open, &tmp, tmp_open)) {
648 dyadic_rational_swap(&tmp, &result.b);
649 result.b_open = tmp_open;
650 }
651
652 lp_dyadic_interval_swap(&result, P);
653 lp_dyadic_interval_destruct(&result);
654 dyadic_rational_destruct(&tmp);
655 }
656 }
657
658 static
lp_value_mul_approx(const lp_value_t * v1,const lp_value_t * v2,lp_value_t * lb,lp_value_t * ub)659 int lp_value_mul_approx(const lp_value_t* v1, const lp_value_t* v2, lp_value_t* lb, lp_value_t* ub) {
660
661 assert(v1 != lb && v1 != ub);
662 assert(v2 != lb && v2 != ub);
663 assert(lb != ub);
664
665 int is_point = 1;
666 lp_integer_t mul_int;
667 lp_rational_t mul_rat;
668 lp_dyadic_rational_t mul_dy;
669 lp_dyadic_interval_t mul_dy_interval;
670
671 if (v1->type == v2->type) {
672 switch (v1->type) {
673 case LP_VALUE_MINUS_INFINITY:
674 case LP_VALUE_PLUS_INFINITY:
675 if (lb) {
676 lp_value_assign_raw(lb, LP_VALUE_PLUS_INFINITY, 0);
677 }
678 if (ub) {
679 lp_value_assign_raw(ub, LP_VALUE_PLUS_INFINITY, 0);
680 }
681 break;
682 case LP_VALUE_INTEGER:
683 lp_integer_construct(&mul_int);
684 lp_integer_mul(lp_Z, &mul_int, &v1->value.z, &v2->value.z);
685 if (lb) {
686 lp_value_assign_raw(lb, LP_VALUE_INTEGER, &mul_int);
687 }
688 if (ub) {
689 lp_value_assign_raw(ub, LP_VALUE_INTEGER, &mul_int);
690 }
691 lp_integer_destruct(&mul_int);
692 break;
693 case LP_VALUE_RATIONAL:
694 lp_rational_construct(&mul_rat);
695 lp_rational_mul(&mul_rat, &v1->value.q, &v2->value.q);
696 if (lb) {
697 lp_value_assign_raw(lb, LP_VALUE_RATIONAL, &mul_rat);
698 }
699 if (ub) {
700 lp_value_assign_raw(ub, LP_VALUE_RATIONAL, &mul_rat);
701 }
702 lp_rational_destruct(&mul_rat);
703 break;
704 case LP_VALUE_DYADIC_RATIONAL:
705 lp_dyadic_rational_construct(&mul_dy);
706 lp_dyadic_rational_mul(&mul_dy, &v1->value.dy_q, &v2->value.dy_q);
707 if (lb) {
708 lp_value_assign_raw(lb, LP_VALUE_DYADIC_RATIONAL, &mul_dy);
709 }
710 if (ub) {
711 lp_value_assign_raw(ub, LP_VALUE_DYADIC_RATIONAL, &mul_dy);
712 }
713 lp_dyadic_rational_destruct(&mul_dy);
714 break;
715 case LP_VALUE_ALGEBRAIC:
716 if (v1->value.a.I.is_point && v2->value.a.I.is_point) {
717 lp_dyadic_rational_construct(&mul_dy);
718 lp_dyadic_rational_mul(&mul_dy, &v1->value.a.I.a, &v2->value.a.I.a);
719 if (lb) {
720 lp_value_assign_raw(lb, LP_VALUE_DYADIC_RATIONAL, &mul_dy);
721 }
722 if (ub) {
723 lp_value_assign_raw(ub, LP_VALUE_DYADIC_RATIONAL, &mul_dy);
724 }
725 lp_dyadic_rational_destruct(&mul_dy);
726 } else {
727 lp_dyadic_interval_construct_zero(&mul_dy_interval);
728 dyadic_interval_mul(&mul_dy_interval, &v1->value.a.I, &v2->value.a.I);
729 if (lb) {
730 lp_value_assign_raw(lb, LP_VALUE_DYADIC_RATIONAL, &mul_dy_interval.a);
731 }
732 if (ub) {
733 lp_value_assign_raw(ub, LP_VALUE_DYADIC_RATIONAL, &mul_dy_interval.b);
734 }
735 is_point = 0;
736 lp_dyadic_interval_destruct(&mul_dy_interval);
737 }
738 break;
739 case LP_VALUE_NONE:
740 assert(0);
741 break;
742 }
743 } else {
744
745 lp_value_t v1_tmp, v2_tmp;
746 const lp_value_t* v1_to_use;
747 const lp_value_t* v2_to_use;
748
749 // Cast to the same type
750 int cast = lp_value_to_same_type(v1, v2, &v1_tmp, &v2_tmp, &v1_to_use, &v2_to_use);
751 if (cast) {
752 // Cast successful, just add as usual
753 lp_value_mul_approx(v1_to_use, v2_to_use, lb, ub);
754 // If a fesh value was used, delete it
755 if (v1_to_use != v1) {
756 lp_value_destruct(&v1_tmp);
757 }
758 if (v2_to_use != v2) {
759 lp_value_destruct(&v2_tmp);
760 }
761 } else {
762 // Cast failed, check the cases
763 int v1_sgn = lp_value_sgn(v1);
764 int v2_sgn = lp_value_sgn(v2);
765
766 if (v1_sgn == 0 || v2_sgn == 0) {
767 if (lb) {
768 lp_value_assign_zero(lb);
769 }
770 if (ub) {
771 lp_value_assign_zero(ub);
772 }
773 return 1;
774 } else {
775 if (lp_value_is_infinity(v1) || lp_value_is_infinity(v2)) {
776 int sgn = v1_sgn*v2_sgn;
777 lp_value_type_t result = LP_VALUE_NONE;
778 if (sgn > 0) {
779 result = LP_VALUE_PLUS_INFINITY;
780 } else if (sgn < 0) {
781 result = LP_VALUE_MINUS_INFINITY;
782 } else {
783 assert(0);
784 }
785 if (lb) {
786 lp_value_assign_raw(lb, result, 0);
787 }
788 if (ub) {
789 lp_value_assign_raw(ub, result, 0);
790 }
791 return 1;
792 }
793 }
794
795 // OK, now deal with algebraic numbers, assume v1 is algebraic
796 if (v2->type == LP_VALUE_ALGEBRAIC) {
797 const lp_value_t* tmp = v1; v1 = v2; v2 = tmp;
798 }
799
800 assert(v1->type == LP_VALUE_ALGEBRAIC);
801 assert(v2->type != LP_VALUE_MINUS_INFINITY);
802 assert(v2->type != LP_VALUE_PLUS_INFINITY);
803 assert(v2->type != LP_VALUE_ALGEBRAIC);
804
805 // v1 is of the form (dy1, dy2), v2 is a proper non-algebraic value
806 // we add them into (dy1 * v2, dy2 * v2). if v2 is negative, we swap the bounds
807
808 lp_value_t v1_lb, v1_ub;
809 const lp_dyadic_rational_t* a = &v1->value.a.I.a;
810 const lp_dyadic_rational_t* b = v1->value.a.I.is_point ? a : &v1->value.a.I.b;
811 lp_value_construct(&v1_lb, LP_VALUE_DYADIC_RATIONAL, a);
812 lp_value_construct(&v1_ub, LP_VALUE_DYADIC_RATIONAL, b);
813 lp_value_mul_approx(&v1_lb, v2, lb, 0);
814 lp_value_mul_approx(&v1_ub, v2, ub, 0);
815 if (v2_sgn < 0) {
816 lp_value_swap(lb, ub);
817 }
818 is_point = 0;
819 lp_value_destruct(&v1_lb);
820 lp_value_destruct(&v1_ub);
821 }
822 }
823
824 return is_point;
825 }
826
lp_interval_mul(lp_interval_t * mul,const lp_interval_t * I1,const lp_interval_t * I2)827 void lp_interval_mul(lp_interval_t* mul, const lp_interval_t* I1, const lp_interval_t* I2) {
828
829 lp_interval_t result;
830 lp_interval_construct_full(&result);
831
832 if (I1->is_point) {
833 if (I2->is_point) {
834 // Just multiply the points
835 result.is_point = lp_value_mul_approx(&I1->a, &I2->a, &result.a, &result.b);
836 if (result.is_point) {
837 lp_value_destruct(&result.b);
838 }
839 result.a_open = result.b_open = !result.is_point;
840 } else {
841 // Depending on the sign of a, we might have to flip
842 int a_sgn = lp_value_sgn(&I1->a);
843 if (a_sgn == 0) {
844 // It's just 0
845 lp_interval_destruct(&result);
846 lp_interval_construct_zero(&result);
847 } else if (a_sgn > 0) {
848 // Regular multiplication
849 int a_is_point = lp_value_mul_approx(&I1->a, &I2->a, &result.a, 0);
850 int b_is_point = lp_value_mul_approx(&I1->a, &I2->b, 0, &result.b);
851 result.is_point = 0;
852 result.a_open = I2->a_open || !a_is_point;
853 result.b_open = I2->b_open || !b_is_point;
854 } else {
855 int a_is_point = lp_value_mul_approx(&I1->a, &I2->b, &result.a, 0);
856 int b_is_point = lp_value_mul_approx(&I1->a, &I2->a, 0, &result.b);
857 result.is_point = 0;
858 result.a_open = I2->b_open || !a_is_point;
859 result.b_open = I2->a_open || !b_is_point;
860 }
861 }
862 } else if (I2->is_point) {
863 return lp_interval_mul(mul, I2, I1);
864 } else {
865
866 assert(!I1->is_point && !I2->is_point);
867
868 //
869 // I1 x I2 = { x*y | x in I1, y in I2 }
870 // = { x*y | I1.a < x < I1.b, I2.a < y < I2.b }
871 //
872
873 lp_value_t tmp_lb, tmp_ub;
874 lp_value_construct_zero(&tmp_lb);
875 lp_value_construct_zero(&tmp_ub);
876
877 int mul_is_point = 0;
878
879 // I1.a x I2.a
880 mul_is_point = lp_value_mul_approx(&I1->a, &I2->a, &result.a, &result.b);
881 result.a_open = result.b_open = I1->a_open || I2->a_open || !mul_is_point;
882
883 // I1.a x I2.b
884 mul_is_point = lp_value_mul_approx(&I1->a, &I2->b, &tmp_lb, &tmp_ub);
885 int tmp_open = I1->a_open || I2->b_open || !mul_is_point;
886 if (lp_interval_endpoint_lt(&tmp_lb, tmp_open, &result.a, result.a_open)) {
887 lp_value_swap(&tmp_lb, &result.a);
888 result.a_open = tmp_open;
889 }
890 if (lp_interval_endpoint_lt(&result.b, result.b_open, &tmp_ub, tmp_open)) {
891 lp_value_swap(&tmp_ub, &result.b);
892 result.b_open = tmp_open;
893 }
894
895 // I1.b x I2.a
896 mul_is_point = lp_value_mul_approx(&I1->b, &I2->a, &tmp_lb, &tmp_ub);
897 tmp_open = I1->b_open || I2->a_open || !mul_is_point;
898 if (lp_interval_endpoint_lt(&tmp_lb, tmp_open, &result.a, result.a_open)) {
899 lp_value_swap(&tmp_lb, &result.a);
900 result.a_open = tmp_open;
901 }
902 if (lp_interval_endpoint_lt(&result.b, result.b_open, &tmp_ub, tmp_open)) {
903 lp_value_swap(&tmp_ub, &result.b);
904 result.b_open = tmp_open;
905 }
906
907 // I1.b x I2.b
908 mul_is_point = lp_value_mul_approx(&I1->b, &I2->b, &tmp_lb, &tmp_ub);
909 tmp_open = I1->b_open || I2->b_open || !mul_is_point;
910 if (lp_interval_endpoint_lt(&tmp_lb, tmp_open, &result.a, result.a_open)) {
911 lp_value_swap(&tmp_lb, &result.a);
912 result.a_open = tmp_open;
913 }
914 if (lp_interval_endpoint_lt(&result.b, result.b_open, &tmp_ub, tmp_open)) {
915 lp_value_swap(&tmp_ub, &result.b);
916 result.b_open = tmp_open;
917 }
918
919 // Final check: if an endpoint is 0, it must have come from another zero
920 // endpoint. If that endpoint is closed, then the resulting endpoint must
921 // be closed too.
922 int sgn_a = lp_value_sgn(&result.a);
923 if (sgn_a == 0) {
924 int c1_a = (lp_value_sgn(&I1->a) == 0) && !I1->a_open;
925 int c1_b = (lp_value_sgn(&I1->b) == 0) && !I1->b_open;
926 int c2_a = (lp_value_sgn(&I2->a) == 0) && !I2->a_open;
927 int c2_b = (lp_value_sgn(&I2->b) == 0) && !I2->b_open;
928 if (c1_a || c1_b || c2_a || c2_b) {
929 result.a_open = 0;
930 }
931 }
932 int sgn_b = lp_value_sgn(&result.b);
933 if (sgn_b == 0) {
934 int c1_a = (lp_value_sgn(&I1->a) == 0) && !I1->a_open;
935 int c1_b = (lp_value_sgn(&I1->b) == 0) && !I1->b_open;
936 int c2_a = (lp_value_sgn(&I2->a) == 0) && !I2->a_open;
937 int c2_b = (lp_value_sgn(&I2->b) == 0) && !I2->b_open;
938 if (c1_a || c1_b || c2_a || c2_b) {
939 result.b_open = 0;
940 }
941 }
942
943 lp_value_destruct(&tmp_lb);
944 lp_value_destruct(&tmp_ub);
945 }
946
947 // Put into result
948 lp_interval_swap(mul, &result);
949 lp_interval_destruct(&result);
950 }
951
952 static
lp_value_pow_approx(const lp_value_t * v,unsigned n,lp_value_t * lb,lp_value_t * ub)953 int lp_value_pow_approx(const lp_value_t* v, unsigned n, lp_value_t* lb, lp_value_t* ub) {
954
955 assert(n > 0);
956 assert(v != lb);
957 assert(v != ub);
958 assert(ub != lb);
959
960 lp_integer_t pow_int;
961 lp_rational_t pow_rat;
962 lp_dyadic_rational_t pow_dy;
963 lp_dyadic_interval_t pow_dy_interval;
964
965 int is_point = 1;
966
967 switch (v->type) {
968 case LP_VALUE_INTEGER:
969 lp_integer_construct(&pow_int);
970 integer_pow(lp_Z, &pow_int, &v->value.z, n);
971 if (lb) {
972 lp_value_assign_raw(lb, LP_VALUE_INTEGER, &pow_int);
973 }
974 if (ub) {
975 lp_value_assign_raw(ub, LP_VALUE_INTEGER, &pow_int);
976 }
977 lp_integer_destruct(&pow_int);
978 break;
979 case LP_VALUE_RATIONAL:
980 lp_rational_construct(&pow_rat);
981 rational_pow(&pow_rat, &v->value.q, n);
982 if (lb) {
983 lp_value_assign_raw(lb, LP_VALUE_RATIONAL, &pow_rat);
984 }
985 if (ub) {
986 lp_value_assign_raw(ub, LP_VALUE_RATIONAL, &pow_rat);
987 }
988 lp_rational_destruct(&pow_rat);
989 break;
990 case LP_VALUE_DYADIC_RATIONAL:
991 lp_dyadic_rational_construct(&pow_dy);
992 dyadic_rational_pow(&pow_dy, &v->value.dy_q, n);
993 if (lb) {
994 lp_value_assign_raw(lb, LP_VALUE_DYADIC_RATIONAL, &pow_dy);
995 }
996 if (ub) {
997 lp_value_assign_raw(ub, LP_VALUE_DYADIC_RATIONAL, &pow_dy);
998 }
999 lp_dyadic_rational_destruct(&pow_dy);
1000 break;
1001 case LP_VALUE_ALGEBRAIC:
1002 if (v->value.a.I.is_point) {
1003 lp_dyadic_rational_construct(&pow_dy);
1004 dyadic_rational_pow(&pow_dy, &v->value.a.I.a, n);
1005 if (lb) {
1006 lp_value_assign_raw(lb, LP_VALUE_DYADIC_RATIONAL, &pow_dy);
1007 }
1008 if (ub) {
1009 lp_value_assign_raw(ub, LP_VALUE_DYADIC_RATIONAL, &pow_dy);
1010 }
1011 lp_dyadic_rational_destruct(&pow_dy);
1012 } else {
1013 lp_dyadic_interval_construct_zero(&pow_dy_interval);
1014 dyadic_interval_pow(&pow_dy_interval, &v->value.a.I, n);
1015 if (lb) {
1016 lp_value_assign_raw(lb, LP_VALUE_DYADIC_RATIONAL, &pow_dy_interval.a);
1017 }
1018 if (ub) {
1019 lp_value_assign_raw(ub, LP_VALUE_DYADIC_RATIONAL, &pow_dy_interval.b);
1020 }
1021 lp_dyadic_interval_destruct(&pow_dy_interval);
1022 is_point = 0;
1023 }
1024 break;
1025 case LP_VALUE_MINUS_INFINITY:
1026 case LP_VALUE_PLUS_INFINITY: {
1027 int sgn = (n % 2) ? lp_value_sgn(v) : 1;
1028 lp_value_type_t type = sgn > 0 ? LP_VALUE_PLUS_INFINITY : LP_VALUE_MINUS_INFINITY;
1029 if (lb) {
1030 lp_value_assign_raw(lb, type, 0);
1031 }
1032 if (ub) {
1033 lp_value_assign_raw(ub, type, 0);
1034 }
1035 break;
1036 }
1037 case LP_VALUE_NONE:
1038 assert(0);
1039 }
1040
1041 return is_point;
1042 }
1043
lp_interval_sgn(const lp_interval_t * I)1044 int lp_interval_sgn(const lp_interval_t* I) {
1045 int a_sgn = lp_value_sgn(&I->a);
1046 if (I->is_point) {
1047 return a_sgn;
1048 }
1049 int b_sgn = lp_value_sgn(&I->b);
1050
1051 if (a_sgn < 0 && b_sgn > 0) {
1052 // Definitively contains 0
1053 return 0;
1054 }
1055 if (a_sgn == 0) {
1056 if (!I->a_open) {
1057 // Left contains zero
1058 return 0;
1059 } else {
1060 return 1;
1061 }
1062 }
1063 if (b_sgn == 0) {
1064 if (!I->b_open) {
1065 // Right contains zero
1066 return 0;
1067 } else {
1068 return -1;
1069 }
1070 }
1071
1072 // Doesn't contain zero
1073 if (a_sgn < 0) {
1074 return -1;
1075 }
1076
1077 assert(b_sgn > 0);
1078 return 1;
1079 }
1080
lp_interval_pow(lp_interval_t * pow,const lp_interval_t * I,unsigned n)1081 void lp_interval_pow(lp_interval_t* pow, const lp_interval_t* I, unsigned n) {
1082
1083 lp_interval_t result;
1084 lp_interval_construct_full(&result);
1085
1086 if (n == 0) {
1087 // I^0 = [1]
1088 lp_value_t one;
1089 lp_value_construct_int(&one, 1);
1090 lp_interval_destruct(&result);
1091 lp_interval_construct_point(&result, &one);
1092 lp_value_destruct(&one);
1093 } else if (I->is_point) {
1094 // Plain power
1095 result.is_point = lp_value_pow_approx(&I->a, n, &result.a, &result.b);
1096 if (result.is_point) {
1097 lp_value_destruct(&result.b);
1098 }
1099 result.a_open = result.b_open = !result.is_point;
1100 } else {
1101 if (n % 2) {
1102 // For odd powers we are monotonic, i.e. [a, b]^n = [a^n, b^n]
1103 int a_point = lp_value_pow_approx(&I->a, n, &result.a, 0);
1104 int b_point = lp_value_pow_approx(&I->b, n, 0, &result.b);
1105 result.a_open = I->a_open || !a_point;
1106 result.b_open = I->b_open || !b_point;
1107 } else {
1108 // Even powers depend on whether 0 is in the interval
1109 int sgn = lp_interval_sgn(I);
1110 if (sgn == 0) {
1111 // P = [0, max(a, b)^n]
1112 int a_point = lp_value_pow_approx(&I->a, n, 0, &result.a);
1113 int b_point = lp_value_pow_approx(&I->b, n, 0, &result.b);
1114 if (lp_interval_endpoint_lt(&result.b, I->b_open, &result.a, I->a_open)) {
1115 lp_value_swap(&result.b, &result.a);
1116 result.b_open = I->a_open || !a_point;
1117 } else {
1118 result.b_open = I->b_open || !b_point;
1119 }
1120 lp_value_assign_zero(&result.a);
1121 result.a_open = 0;
1122 } else if (sgn > 0) {
1123 // P = I^n
1124 int a_point = lp_value_pow_approx(&I->a, n, &result.a, 0);
1125 int b_point = lp_value_pow_approx(&I->b, n, 0, &result.b);
1126 result.a_open = I->a_open || !a_point;
1127 result.b_open = I->b_open || !b_point;
1128 } else {
1129 // P = I^n, but swappeed
1130 int a_point = lp_value_pow_approx(&I->a, n, 0, &result.b);
1131 int b_point = lp_value_pow_approx(&I->b, n, &result.a, 0);
1132 result.a_open = I->b_open || !b_point;
1133 result.b_open = I->a_open || !a_point;
1134 }
1135 }
1136 }
1137
1138 lp_interval_swap(pow, &result);
1139 lp_interval_destruct(&result);
1140 }
1141
rational_interval_pow(lp_rational_interval_t * P,const lp_rational_interval_t * I,unsigned n)1142 void rational_interval_pow(lp_rational_interval_t* P, const lp_rational_interval_t* I, unsigned n) {
1143 if (n == 0) {
1144 // I^0 = [1]
1145 if (!P->is_point) {
1146 P->is_point = 1;
1147 rational_destruct(&P->b);
1148 }
1149 rational_assign_int(&P->a, 1, 1);
1150 P->a_open = 0;
1151 P->b_open = 0;
1152 } else if (I->is_point) {
1153 // Plain power
1154 if (!P->is_point) {
1155 rational_destruct(&P->b);
1156 P->is_point = 1;
1157 P->a_open = P->b_open = 0;
1158 }
1159 rational_pow(&P->a, &I->a, n);
1160 } else {
1161 if (P->is_point) {
1162 P->is_point = 0;
1163 rational_construct(&P->b);
1164 }
1165 if (n % 2) {
1166 // For odd powers we are monotonic, i.e. [a, b]^n = [a^n, b^n]
1167 P->a_open = I->a_open;
1168 P->b_open = I->b_open;
1169 rational_pow(&P->a, &I->a, n);
1170 rational_pow(&P->b, &I->b, n);
1171 } else {
1172 // Even powers depend on whether 0 is in the interval
1173 int sgn = lp_rational_interval_sgn(I);
1174 rational_pow(&P->a, &I->a, n);
1175 rational_pow(&P->b, &I->b, n);
1176 if (sgn == 0) {
1177 // P = [0, max(a, b)^n]
1178 if (rational_interval_endpoint_lt(&P->b, I->b_open, &P->a, I->a_open)) {
1179 rational_swap(&P->b, &P->a);
1180 P->b_open = I->a_open;
1181 } else {
1182 P->b_open = I->b_open;
1183 }
1184 rational_assign_int(&P->a, 0, 1);
1185 P->a_open = 0;
1186 } else if (sgn > 0) {
1187 // P = I^n
1188 P->a_open = I->a_open;
1189 P->b_open = I->b_open;
1190 } else {
1191 // negative turns positive, so we flip
1192 rational_swap(&P->a, &P->b);
1193 P->a_open = I->b_open;
1194 P->b_open = I->a_open;
1195 }
1196 }
1197 }
1198 }
1199
dyadic_interval_pow(lp_dyadic_interval_t * P,const lp_dyadic_interval_t * I,unsigned n)1200 void dyadic_interval_pow(lp_dyadic_interval_t* P, const lp_dyadic_interval_t* I, unsigned n) {
1201 if (n == 0) {
1202 // I^0 = [1]
1203 if (!P->is_point) {
1204 P->is_point = 1;
1205 dyadic_rational_destruct(&P->b);
1206 }
1207 dyadic_rational_assign_int(&P->a, 1, 1);
1208 P->a_open = 0;
1209 P->b_open = 0;
1210 } else if (I->is_point) {
1211 // Plain power
1212 if (!P->is_point) {
1213 dyadic_rational_destruct(&P->b);
1214 P->is_point = 1;
1215 P->a_open = P->b_open = 0;
1216 }
1217 dyadic_rational_pow(&P->a, &I->a, n);
1218 } else {
1219 if (P->is_point) {
1220 P->is_point = 0;
1221 dyadic_rational_construct(&P->b);
1222 }
1223 if (n % 2) {
1224 // For odd powers we are monotonic, i.e. [a, b]^n = [a^n, b^n]
1225 P->a_open = I->a_open;
1226 P->b_open = I->b_open;
1227 dyadic_rational_pow(&P->a, &I->a, n);
1228 dyadic_rational_pow(&P->b, &I->b, n);
1229 } else {
1230 // Even powers depend on whether 0 is in the interval
1231 int sgn = lp_dyadic_interval_sgn(I);
1232 dyadic_rational_pow(&P->a, &I->a, n);
1233 dyadic_rational_pow(&P->b, &I->b, n);
1234 if (sgn == 0) {
1235 // P = [0, max(a, b)^n]
1236 if (dyadic_interval_endpoint_lt(&P->b, I->b_open, &P->a, I->a_open)) {
1237 dyadic_rational_swap(&P->b, &P->a);
1238 P->b_open = I->a_open;
1239 } else {
1240 P->b_open = I->b_open;
1241 }
1242 dyadic_rational_assign_int(&P->a, 0, 1);
1243 P->a_open = 0;
1244 } else if (sgn > 0) {
1245 // P = I^n
1246 P->a_open = I->a_open;
1247 P->b_open = I->b_open;
1248 } else {
1249 // negative turns positive, so we flip
1250 dyadic_rational_swap(&P->a, &P->b);
1251 P->a_open = I->b_open;
1252 P->b_open = I->a_open;
1253 }
1254 }
1255 }
1256 }
1257