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