1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* File....: term2.c */
4 /* Name....: Term Functions */
5 /* Author..: Thorsten Koch */
6 /* Copyright by Author, All rights reserved */
7 /* */
8 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
9 /*
10 * Copyright (C) 2001-2018 by Thorsten Koch <koch@zib.de>
11 *
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 3
15 * of the License, or (at your option) any later version.
16 *
17 * This program is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public License
23 * along with this program; if not, write to the Free Software
24 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
25 */
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <string.h>
30 #include <assert.h>
31
32 /* #define TRACE 1 */
33
34 #include <stdbool.h>
35 #include "zimpl/mshell.h"
36 #include "zimpl/ratlptypes.h"
37 #include "zimpl/numb.h"
38 #include "zimpl/elem.h"
39 #include "zimpl/tuple.h"
40 #include "zimpl/mme.h"
41 #include "zimpl/bound.h"
42 #include "zimpl/set.h"
43 #include "zimpl/entry.h"
44 #include "zimpl/mono.h"
45 #include "zimpl/hash.h"
46 #include "zimpl/term.h"
47 #include "zimpl/stmt.h"
48 #include "zimpl/prog.h"
49 #include "zimpl/xlpglue.h"
50
51
52 #define TERM_EXTEND_SIZE 16
53 #define TERM_SID 0x5465726d
54
55 struct term
56 {
57 SID
58 Numb* constant;
59 int size;
60 int used;
61 Mono** elem;
62 };
63
term_new(int size)64 Term* term_new(int size)
65 {
66 Term* term = calloc(1, sizeof(*term));
67
68 Trace("term_new");
69
70 assert(term != NULL);
71 assert(size > 0);
72
73 term->constant = numb_new_integer(0);
74 term->size = size;
75 term->used = 0;
76 term->elem = calloc((size_t)term->size, sizeof(*term->elem));
77
78 SID_set(term, TERM_SID);
79 assert(term_is_valid(term));
80
81 return term;
82 }
83
term_add_elem(Term * term,const Entry * entry,const Numb * coeff,MFun mfun)84 void term_add_elem(Term* term, const Entry* entry, const Numb* coeff, MFun mfun)
85 {
86 Trace("term_add_elem");
87
88 assert(term_is_valid(term));
89 assert(entry_is_valid(entry));
90 assert(!numb_equal(coeff, numb_zero()));
91 assert(term->used <= term->size);
92
93 if (term->used == term->size)
94 {
95 term->size += TERM_EXTEND_SIZE;
96 term->elem = realloc(
97 term->elem, (size_t)term->size * sizeof(*term->elem));
98
99 assert(term->elem != NULL);
100 }
101 assert(term->used < term->size);
102
103 term->elem[term->used] = mono_new(coeff, entry, mfun);
104 term->used++;
105
106 assert(term_is_valid(term));
107 }
108
109 #if 0 /* not used */
110 void term_mul_elem(Term* term, const Entry* entry, const Numb* coeff)
111 {
112 int i;
113
114 Trace("term_mul_elem");
115
116 assert(term_is_valid(term));
117 assert(entry_is_valid(entry));
118 assert(!numb_equal(coeff, numb_zero()));
119 assert(term->used <= term->size);
120
121 for(i = 0; i < term->used; i++)
122 {
123 mono_mul_entry(term->elem[i], entry);
124 mono_mul_coeff(term->elem[i], coeff);
125 }
126 assert(term_is_valid(term));
127 }
128 #endif
129
term_free(Term * term)130 void term_free(Term* term)
131 {
132 int i;
133
134 Trace("term_free");
135
136 assert(term_is_valid(term));
137
138 SID_del(term);
139
140 for(i = 0; i < term->used; i++)
141 mono_free(term->elem[i]);
142
143 numb_free(term->constant);
144
145 free(term->elem);
146 free(term);
147 }
148
term_is_valid(const Term * term)149 bool term_is_valid(const Term* term)
150 {
151 int i;
152
153 if (term == NULL || !SID_ok(term, TERM_SID) || term->used > term->size)
154 abort();
155
156 for(i = 0; i < term->used; i++)
157 if (numb_equal(mono_get_coeff(term->elem[i]), numb_zero()))
158 abort();
159
160 return true;
161 }
162
term_copy(const Term * term)163 Term* term_copy(const Term* term)
164 {
165 Term* tnew = term_new(term->used + TERM_EXTEND_SIZE);
166 int i;
167
168 Trace("term_copy");
169
170 assert(term_is_valid(term));
171 assert(term_is_valid(tnew));
172
173 for(i = 0; i < term->used; i++)
174 tnew->elem[i] = mono_copy(term->elem[i]);
175
176 tnew->used = term->used;
177 numb_set(tnew->constant, term->constant);
178
179 assert(term_is_valid(tnew));
180
181 return tnew;
182 }
183
term_append_term(Term * term,const Term * term_b)184 void term_append_term(
185 Term* term,
186 const Term* term_b)
187 {
188 int i;
189
190 /* ??? test auf gleiche monome fehlt!!! */
191
192 Trace("term_append_term");
193
194 assert(term_is_valid(term));
195 assert(term_is_valid(term_b));
196
197 if (term->used + term_b->used >= term->size)
198 {
199 term->size = term->used + term_b->used;
200 term->elem = realloc(
201 term->elem, (size_t)(term->size + TERM_EXTEND_SIZE)* sizeof(*term->elem));
202
203 assert(term->elem != NULL);
204 }
205 assert(term->used + term_b->used <= term->size);
206
207 if (!numb_equal(term_b->constant, numb_zero()))
208 numb_add(term->constant, term_b->constant);
209
210 for(i = 0; i < term_b->used; i++)
211 {
212 term->elem[term->used] = mono_copy(term_b->elem[i]);
213 term->used++;
214 }
215 assert(term_is_valid(term));
216 }
217
term_mul_term(const Term * term_a,const Term * term_b)218 Term* term_mul_term(const Term* term_a, const Term* term_b)
219 {
220 Term* term;
221 Term* term_simplified;
222 int i;
223 int k;
224
225 Trace("term_mul_term");
226
227 assert(term_is_valid(term_a));
228 assert(term_is_valid(term_b));
229
230 term = term_new((term_a->used + 1) * (term_b->used + 1)); /* +1 for constant */
231
232 for(i = 0; i < term_a->used; i++)
233 {
234 for(k = 0; k < term_b->used; k++)
235 {
236 assert(term->used < term->size);
237
238 term->elem[term->used] = mono_mul(term_a->elem[i], term_b->elem[k]);
239 term->used++;
240 }
241 }
242 if (!numb_equal(term_b->constant, numb_zero()))
243 {
244 for(i = 0; i < term_a->used; i++)
245 {
246 assert(term->used < term->size);
247
248 term->elem[term->used] = mono_copy(term_a->elem[i]);
249 mono_mul_coeff(term->elem[term->used], term_b->constant);
250 term->used++;
251 }
252 }
253 if (!numb_equal(term_a->constant, numb_zero()))
254 {
255 for(i = 0; i < term_b->used; i++)
256 {
257 assert(term->used < term->size);
258
259 term->elem[term->used] = mono_copy(term_b->elem[i]);
260 mono_mul_coeff(term->elem[term->used], term_a->constant);
261 term->used++;
262 }
263 }
264 numb_free(term->constant);
265 term->constant = numb_new_mul(term_a->constant, term_b->constant);
266
267 assert(term_is_valid(term));
268
269 term_simplified = term_simplify(term);
270
271 term_free(term);
272
273 return term_simplified;
274 }
275
term_add_term(const Term * term_a,const Term * term_b)276 Term* term_add_term(const Term* term_a, const Term* term_b)
277 {
278 Term* term;
279 int i;
280
281 Trace("term_add_term");
282
283 assert(term_is_valid(term_a));
284 assert(term_is_valid(term_b));
285
286 term = term_new(term_a->used + term_b->used + TERM_EXTEND_SIZE);
287 term->used = term_a->used + term_b->used;
288
289 numb_set(term->constant, term_a->constant);
290 numb_add(term->constant, term_b->constant);
291
292 assert(term->size >= term->used);
293
294 for(i = 0; i < term_a->used; i++)
295 term->elem[i] = mono_copy(term_a->elem[i]);
296
297 for(i = 0; i < term_b->used; i++)
298 term->elem[i + term_a->used] = mono_copy(term_b->elem[i]);
299
300 assert(term_is_valid(term));
301
302 return term;
303 }
304
term_sub_term(const Term * term_a,const Term * term_b)305 Term* term_sub_term(const Term* term_a, const Term* term_b)
306 {
307 Term* term;
308 int i;
309
310 Trace("term_sub_term");
311
312 assert(term_is_valid(term_a));
313 assert(term_is_valid(term_b));
314
315 term = term_new(term_a->used + term_b->used + TERM_EXTEND_SIZE);
316 term->used = term_a->used + term_b->used;
317
318 numb_set(term->constant, term_a->constant);
319 numb_sub(term->constant, term_b->constant);
320
321 assert(term->size >= term->used);
322
323 for(i = 0; i < term_a->used; i++)
324 term->elem[i] = mono_copy(term_a->elem[i]);
325
326 for(i = 0; i < term_b->used; i++)
327 {
328 term->elem[i + term_a->used] = mono_copy(term_b->elem[i]);
329 mono_neg(term->elem[i + term_a->used]);
330 }
331 assert(term_is_valid(term));
332
333 return term;
334 }
335
336 /** Combines monoms in the term where possible
337 */
338 /* ??? TODO:reimplement with a hash list */
339 #if 1
term_simplify(const Term * term_org)340 Term* term_simplify(const Term* term_org)
341 {
342 Term* term;
343 int i;
344 Hash* hash;
345
346 assert(term_is_valid(term_org));
347
348 term = term_new(term_org->used);
349 hash = hash_new(HASH_MONO, term_org->used);
350
351 numb_set(term->constant, term_org->constant);
352
353 for(i = 0; i < term_org->used; i++)
354 {
355 const Mono* mono = hash_lookup_mono(hash, term_org->elem[i]);
356
357 if (mono == NULL)
358 {
359 term->elem[term->used] = mono_copy(term_org->elem[i]);
360
361 hash_add_mono(hash, term->elem[term->used]);
362
363 term->used++;
364 }
365 else
366 {
367 assert(mono_equal(mono, term_org->elem[i]));
368
369 mono_add_coeff(mono, mono_get_coeff(term_org->elem[i]));
370 }
371 }
372 hash_free(hash);
373
374 /* Check if there are any cancellations
375 */
376 for(i = 0; i < term->used; i++)
377 {
378 if (numb_equal(mono_get_coeff(term->elem[i]), numb_zero()))
379 {
380 mono_free(term->elem[i]);
381
382 if (term->used > 0)
383 {
384 term->used--;
385 term->elem[i] = term->elem[term->used];
386 }
387 }
388 }
389 assert(term_is_valid(term));
390
391 return term;
392 }
393 #else /* Old and quadratic */
term_simplify(const Term * term_org)394 Term* term_simplify(const Term* term_org)
395 {
396 Term* term;
397 Bool* done;
398 int i;
399
400 assert(term_is_valid(term_org));
401
402 term = term_new(term_org->used);
403 done = calloc(term_org->used, sizeof(*done));
404
405 assert(done != NULL);
406
407 numb_set(term->constant, term_org->constant);
408
409 for(i = 0; i < term_org->used; i++)
410 {
411 int k;
412
413 if (done[i])
414 continue;
415
416 term->elem[term->used] = mono_copy(term_org->elem[i]);
417
418 for(k = i + 1; k < term_org->used; k++)
419 {
420 if (!done[k] && mono_equal(term->elem[term->used], term_org->elem[k]))
421 {
422 mono_add_coeff(term->elem[term->used], mono_get_coeff(term_org->elem[k]));
423 done[k] = true;
424 }
425 }
426 if (!numb_equal(mono_get_coeff(term->elem[term->used]), numb_zero()))
427 term->used++;
428 else
429 mono_free(term->elem[term->used]);
430 }
431 free(done);
432
433 assert(term_is_valid(term));
434
435 return term;
436 }
437 #endif
438
439 /*lint -e{818} supress "Pointer parameter 'term' could be declared as pointing to const" */
term_add_constant(Term * term,const Numb * value)440 void term_add_constant(
441 Term* term,
442 const Numb* value)
443 {
444 Trace("term_add_constant");
445
446 assert(term_is_valid(term));
447
448 numb_add(term->constant, value);
449
450 assert(term_is_valid(term));
451 }
452
453 /*lint -e{818} supress "Pointer parameter 'term' could be declared as pointing to const" */
term_sub_constant(Term * term,const Numb * value)454 void term_sub_constant(Term* term, const Numb* value)
455 {
456 Trace("term_sub_constant");
457
458 assert(term_is_valid(term));
459
460 numb_sub(term->constant, value);
461
462 assert(term_is_valid(term));
463 }
464
term_mul_coeff(Term * term,const Numb * value)465 void term_mul_coeff(Term* term, const Numb* value)
466 {
467 int i;
468
469 Trace("term_mul_coeff");
470
471 assert(term_is_valid(term));
472
473 numb_mul(term->constant, value);
474
475 if (numb_equal(value, numb_zero()))
476 {
477 for(i = 0; i < term->used; i++)
478 mono_free(term->elem[i]);
479
480 term->used = 0;
481 }
482 else
483 {
484 for(i = 0; i < term->used; i++)
485 mono_mul_coeff(term->elem[i], value);
486 }
487 assert(term_is_valid(term));
488 }
489
term_get_constant(const Term * term)490 const Numb* term_get_constant(const Term* term)
491 {
492 assert(term_is_valid(term));
493
494 return term->constant;
495 }
496
497 #if 0 /* not used */
498 void term_negate(Term* term)
499 {
500 Trace("term_negate");
501
502 assert(term_is_valid(term));
503
504 term_mul_coeff(term, numb_minusone());
505 }
506 #endif
507
term_to_objective(const Term * term)508 void term_to_objective(const Term* term)
509 {
510 Var* var;
511 int i;
512
513 Trace("term_to_objective");
514
515 assert(term_is_valid(term));
516
517 /* If we have a constant in the objective, generate an artificial
518 * variable @OOffset fixed to the value. This works allways.
519 * This is needed because there is no general way for example in
520 * MPS format to specify an objective value offset.
521 */
522 if (!numb_equal(term->constant, numb_zero()))
523 {
524 const char* format = "%sObjOffset";
525
526 Bound* lower = bound_new(BOUND_VALUE, numb_one());
527 Bound* upper = bound_new(BOUND_VALUE, numb_one());
528 char* vname = malloc(strlen(SYMBOL_NAME_INTERNAL) + strlen(format) + 1);
529
530 sprintf(vname, format, SYMBOL_NAME_INTERNAL);
531 var = xlp_addvar(prog_get_lp(), vname, VAR_CON, lower, upper, numb_zero(), numb_zero());
532 xlp_addtocost(prog_get_lp(), var, term->constant);
533
534 free(vname);
535 bound_free(upper);
536 bound_free(lower);
537 }
538
539 for(i = 0; i < term->used; i++)
540 {
541 const Numb* coeff = mono_get_coeff(term->elem[i]);
542
543 assert(!numb_equal(coeff, numb_zero()));
544 assert(mono_is_linear(term->elem[i]));
545
546 var = mono_get_var(term->elem[i], 0);
547
548 xlp_addtocost(prog_get_lp(), var, coeff);
549 }
550 }
551
term_get_elements(const Term * term)552 int term_get_elements(const Term* term)
553 {
554 assert(term_is_valid(term));
555
556 return term->used;
557 }
558
term_get_element(const Term * term,int i)559 Mono* term_get_element(const Term* term, int i)
560 {
561 assert(term_is_valid(term));
562 assert(i >= 0);
563 assert(i < term->used);
564
565 return term->elem[i];
566 }
567
term_get_lower_bound(const Term * term)568 Bound* term_get_lower_bound(const Term* term)
569 {
570 Bound* bound;
571 Numb* lower;
572 Numb* value;
573 int i;
574
575 lower = numb_new_integer(0);
576
577 for(i = 0; i < term->used; i++)
578 {
579 const Numb* coeff = mono_get_coeff(term->elem[i]);
580 int sign = numb_get_sgn(coeff);
581
582 assert(sign != 0);
583
584 bound = (sign > 0)
585 ? xlp_getlower(prog_get_lp(), mono_get_var(term->elem[i], 0))
586 : xlp_getupper(prog_get_lp(), mono_get_var(term->elem[i], 0));
587
588 if (bound_get_type(bound) != BOUND_VALUE)
589 goto finish;
590
591 value = numb_new_mul(bound_get_value(bound), coeff);
592
593 numb_add(lower, value);
594
595 bound_free(bound);
596 numb_free(value);
597 }
598 numb_add(lower, term->constant);
599
600 bound = bound_new(BOUND_VALUE, lower);
601
602 finish:
603 numb_free(lower);
604
605 return bound;
606 }
607
term_get_upper_bound(const Term * term)608 Bound* term_get_upper_bound(const Term* term)
609 {
610 Bound* bound;
611 Numb* upper;
612 Numb* value;
613 int i;
614
615 upper = numb_new_integer(0);
616
617 for(i = 0; i < term->used; i++)
618 {
619 const Numb* coeff = mono_get_coeff(term->elem[i]);
620 int sign = numb_get_sgn(coeff);
621
622 assert(sign != 0);
623
624 bound = (sign > 0)
625 ? xlp_getupper(prog_get_lp(), mono_get_var(term->elem[i], 0))
626 : xlp_getlower(prog_get_lp(), mono_get_var(term->elem[i], 0));
627
628 if (bound_get_type(bound) != BOUND_VALUE)
629 goto finish;
630
631 value = numb_new_mul(bound_get_value(bound), coeff);
632
633 numb_add(upper, value);
634
635 bound_free(bound);
636 numb_free(value);
637 }
638 numb_add(upper, term->constant);
639
640 bound = bound_new(BOUND_VALUE, upper);
641
642 finish:
643 numb_free(upper);
644
645 return bound;
646 }
647
term_get_degree(const Term * term)648 int term_get_degree(const Term* term)
649 {
650 int degree = 0;
651 int i;
652
653 for(i = 0; i < term->used; i++)
654 if (degree < mono_get_degree(term->elem[i]))
655 degree = mono_get_degree(term->elem[i]);
656
657 return degree;
658 }
659
term_is_linear(const Term * term)660 bool term_is_linear(const Term* term)
661 {
662 int i;
663
664 for(i = 0; i < term->used; i++)
665 if (!mono_is_linear(term->elem[i]))
666 return false;
667
668 return true;
669 }
670
term_is_polynomial(const Term * term)671 bool term_is_polynomial(const Term* term)
672 {
673 int i;
674
675 for(i = 0; i < term->used; i++)
676 if ((mono_get_function(term->elem[i]) != MFUN_NONE)
677 && (mono_get_function(term->elem[i]) != MFUN_TRUE)
678 && (mono_get_function(term->elem[i]) != MFUN_FALSE))
679 return false;
680
681 return true;
682 }
683
684
term_make_conditional(const Term * ind_term,const Term * cond_term,bool is_true)685 Term* term_make_conditional(const Term* ind_term, const Term* cond_term, bool is_true)
686 {
687 Term* term;
688
689 assert(term_is_valid(ind_term));
690 assert(term_is_valid(cond_term));
691
692 assert(term_get_elements(ind_term) == 1);
693 assert(term_get_degree(ind_term) == 1);
694 assert(numb_equal(mono_get_coeff(term_get_element(ind_term, 0)), numb_one()));
695
696 term = term_copy(ind_term);
697
698 mono_set_function(term_get_element(term, 0), is_true ? MFUN_TRUE : MFUN_FALSE);
699
700 term_append_term(term, cond_term);
701
702 return term;
703 }
704
term_is_all_integer(const Term * term)705 bool term_is_all_integer(const Term* term)
706 {
707 int i;
708
709 for(i = 0; i < term->used; i++)
710 {
711 VarClass vc = xlp_getclass(prog_get_lp(), mono_get_var(term->elem[i], 0));
712
713 if (vc != VAR_INT && vc != VAR_IMP)
714 return false;
715 }
716 return true;
717 }
718
719 #ifndef NDEBUG
term_print(FILE * fp,const Term * term,bool print_symbol_index)720 void term_print(FILE* fp, const Term* term, bool print_symbol_index)
721 {
722 int i;
723
724 assert(term_is_valid(term));
725
726 for(i = 0; i < term->used; i++)
727 mono_print(fp, term->elem[i], print_symbol_index);
728
729 if (!numb_equal(term->constant, numb_zero()))
730 {
731 if (numb_cmp(term->constant, numb_zero()) >= 0)
732 fprintf(fp, " + %.16g ", numb_todbl(term->constant));
733 else
734 fprintf(fp, " - %.16g ", -numb_todbl(term->constant));
735 }
736 }
737 #endif /* !NDEBUG */
738
739
740
741
742
743
744
745