1 /*
2 This file is part of PolyLib.
3
4 PolyLib is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 PolyLib is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with PolyLib. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 /* vector.c
19 COPYRIGHT
20 Both this software and its documentation are
21
22 Copyrighted 1993 by IRISA /Universite de Rennes I - France,
23 Copyright 1995,1996 by BYU, Provo, Utah
24 all rights reserved.
25
26 Permission is granted to copy, use, and distribute
27 for any commercial or noncommercial purpose under the terms
28 of the GNU General Public license, version 2, June 1991
29 (see file : LICENSING).
30 */
31
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include <polylib/polylib.h>
35
36 #ifdef MAC_OS
37 #define abs __abs
38 #endif
39
40 /*
41 * Compute n!
42 */
43 void Factorial(int n, Value *fact) {
44
45 int i;
46 Value tmp;
47
48 value_init(tmp);
49
50 value_set_si(*fact,1);
51 for (i=1;i<=n;i++) {
52 value_set_si(tmp,i);
53 value_multiply(*fact,*fact,tmp);
54 }
55 value_clear(tmp);
56 } /* Factorial */
57
58 /*
59 * Compute n!/(p!(n-p)!)
60 */
61 void Binomial(int n, int p, Value *result) {
PDomainIntersection(Polyhedron * Pol1,Polyhedron * Pol2,unsigned NbMaxRays)62
63 int i;
64 Value tmp;
65 Value f;
66
67 value_init(tmp);value_init(f);
68
69 if (n-p<p)
70 p=n-p;
71 if (p!=0) {
72 value_set_si(*result,(n-p+1));
73 for (i=n-p+2;i<=n;i++) {
74 value_set_si(tmp,i);
75 value_multiply(*result,*result,tmp);
76 }
77 Factorial(p,&f);
78 value_division(*result,*result,f);
79 }
80 else
81 value_set_si(*result,1);
82 value_clear(f);value_clear(tmp);
83 } /* Binomial */
84
85 /*
86 * Return the number of ways to choose 'b' items from 'a' items, that is,
87 * return a!/(b!(a-b)!)
88 */
89 void CNP(int a,int b, Value *result) {
90
91 int i;
92 Value tmp;
93 value_init(tmp);
94
95 value_set_si(*result,1);
96
97 /* If number of items is less than the number to be choosen, return 1 */
98 if(a <= b) {
99 value_clear(tmp);
100 return;
101 }
102 for(i=a;i>b;--i) {
PDomainDifference(Polyhedron * Pol1,Polyhedron * Pol2,unsigned NbMaxRays)103 value_set_si(tmp,i);
104 value_multiply(*result,*result,tmp);
105 }
106 for(i=1;i<=(a-b);++i) {
107 value_set_si(tmp,i);
108 value_division(*result,*result,tmp);
109 }
110 value_clear(tmp);
111 } /* CNP */
112
113 /*
114 * Compute GCD of 'a' and 'b'
115 */
116 void Gcd(Value a,Value b,Value *result) {
117
118 Value acopy, bcopy;
119
120 value_init(acopy);
121 value_init(bcopy);
122 value_assign(acopy,a);
123 value_assign(bcopy,b);
124 while(value_notzero_p(acopy)) {
125 value_modulus(*result,bcopy,acopy);
126 value_assign(bcopy,acopy);
127 value_assign(acopy,*result);
128 }
129 value_absolute(*result,bcopy);
130 value_clear(acopy);
131 value_clear(bcopy);
132 } /* Gcd */
133
134 /*
135 * Return the smallest component index in 'p' whose value is non-zero
136 */
137 int First_Non_Zero(Value *p,unsigned length) {
138
139 Value *cp;
140 int i;
141
142 cp = p;
143 for (i=0;i<length;i++) {
144 if (value_notzero_p(*cp))
145 break;
146 cp++;
147 }
148 return((i==length) ? -1 : i );
149 } /* First_Non_Zero */
150
TestRank(Matrix * Mat)151 /*
152 * Allocate memory space for Vector
153 */
154 Vector *Vector_Alloc(unsigned length) {
155
156 int i;
157 Vector *vector;
158
159 vector = (Vector *)malloc(sizeof(Vector));
160 if (!vector) {
161 errormsg1("Vector_Alloc", "outofmem", "out of memory space");
162 return 0;
163 }
164 vector->Size=length;
165 vector->p=(Value *)malloc(length * sizeof(Value));
166 if (!vector->p) {
167 errormsg1("Vector_Alloc", "outofmem", "out of memory space");
168 free(vector);
169 return 0;
170 }
171 for(i=0;i<length;i++)
172 value_init(vector->p[i]);
173 return vector;
174 } /* Vector_Alloc */
175
176 /*
177 * Free the memory space occupied by Vector
178 */
179 void Vector_Free(Vector *vector) {
180
181 int i;
182
183 if (!vector) return;
184 for(i=0;i<vector->Size;i++)
185 value_clear(vector->p[i]);
186 free(vector->p);
187 free(vector);
188 } /* Vector_Free */
189
190 /*
191 * Print the contents of a Vector
192 */
193 void Vector_Print(FILE *Dst, const char *Format, Vector *vector)
194 {
195 int i;
196 Value *p;
197 unsigned length;
198
199 fprintf(Dst, "%d\n", length=vector->Size);
200 p = vector->p;
201 for (i=0;i<length;i++) {
202 if (Format) {
203 value_print(Dst,Format,*p++);
204 }
205 else {
206 value_print(Dst,P_VALUE_FMT,*p++);
207 }
208 }
209 fprintf(Dst, "\n");
210 } /* Vector_Print */
211
212 /*
213 * Read the contents of a Vector
214 */
215 Vector *Vector_Read() {
216
217 Vector *vector;
218 unsigned length;
219 int i;
220 char str[1024];
221 Value *p;
222
223 scanf("%d", &length);
224 vector = Vector_Alloc(length);
225 if (!vector) {
226 errormsg1("Vector_Read", "outofmem", "out of memory space");
227 return 0;
228 }
SMAlloc(int rows,int cols)229 p = vector->p;
230 for (i=0;i<length;i++) {
231 scanf("%s",str);
232 value_read(*(p++),str);
233 }
234 return vector;
235 } /* Vector_Read */
236
237 /*
238 * Assign 'n' to each component of Vector 'p'
239 */
240 void Vector_Set(Value *p,int n,unsigned length) {
241
242 Value *cp;
243 int i;
244
245 cp = p;
246 for (i=0;i<length;i++) {
247 value_set_si(*cp,n);
248 cp++;
249 }
250 return;
251 } /* Vector_Set */
SMPrint(SatMatrix * matrix)252
253 /*
254 * Exchange the components of the vectors 'p1' and 'p2' of length 'length'
255 */
256 void Vector_Exchange(Value *p1, Value *p2, unsigned length) {
257
258 int i;
259
260 for(i=0;i<length;i++) {
261 value_swap(p1[i],p2[i]);
262 }
263 return;
264 }
265
266 /*
267 * Copy Vector 'p1' to Vector 'p2'
SMFree(SatMatrix * matrix)268 */
269 void Vector_Copy(Value *p1,Value *p2,unsigned length) {
270
271 int i;
272 Value *cp1, *cp2;
273
274 cp1 = p1;
275 cp2 = p2;
276
277 for(i=0;i<length;i++)
278 value_assign(*cp2++,*cp1++);
279
280 return;
281 }
282
283 /*
284 * Add two vectors 'p1' and 'p2' and store the result in 'p3'
285 */
286 void Vector_Add(Value *p1,Value *p2,Value *p3,unsigned length) {
287
288 Value *cp1, *cp2, *cp3;
289 int i;
290
291 cp1=p1;
292 cp2=p2;
293 cp3=p3;
294 for (i=0;i<length;i++) {
295
296 /* *cp3++ = *cp1++ + *cp2++ */
297 value_addto(*cp3,*cp1,*cp2);
298 cp1++; cp2++; cp3++;
299 }
300 } /* Vector_Add */
301
302 /*
303 * Subtract two vectors 'p1' and 'p2' and store the result in 'p3'
304 */
305 void Vector_Sub(Value *p1,Value *p2,Value *p3,unsigned length) {
306
307 Value *cp1, *cp2, *cp3;
308 int i;
309
310 cp1=p1;
311 cp2=p2;
312 cp3=p3;
313 for (i=0;i<length;i++) {
Add_CEqualities(Polyhedron * D)314
315 /* *cp3++= *cp1++ - *cp2++ */
316 value_subtract(*cp3,*cp1,*cp2);
317 cp1++; cp2++; cp3++;
318 }
319 } /* Vector_Sub */
320
321 /*
322 * Compute Bitwise OR of Vectors 'p1' and 'p2' and store in 'p3'
323 */
324 void Vector_Or(Value *p1,Value *p2,Value *p3,unsigned length) {
325
326 Value *cp1, *cp2, *cp3;
327 int i;
328
329 cp1=p1;
330 cp2=p2;
331 cp3=p3;
332 for (i=0;i<length;i++) {
333
334 /* *cp3++=*cp1++ | *cp2++ */
335 value_orto(*cp3,*cp1,*cp2);
336 cp1++; cp2++; cp3++;
337 }
338 } /* Vector_Or */
339
340 /*
int_array2bit_vector(unsigned int * array,int n)341 * Scale Vector 'p1' lambda times and store in 'p2'
342 */
343 void Vector_Scale(Value *p1,Value *p2,Value lambda,unsigned length) {
344
345 Value *cp1, *cp2;
346 int i;
347
348 cp1=p1;
349 cp2=p2;
350 for (i=0;i<length;i++) {
351
352 /* *cp2++=*cp1++ * lambda */
353 value_multiply(*cp2,*cp1,lambda);
354 cp1++; cp2++;
355 }
356 } /* Vector_Scale */
357
358 /*
359 * Antiscale Vector 'p1' by lambda and store in 'p2'
360 * Assumes all elements of 'p1' are divisble by lambda.
361 */
362 void Vector_AntiScale(Value *p1, Value *p2, Value lambda, unsigned length)
traite_m_face(Polyhedron * D,unsigned int * mf,unsigned int * egalite)363 {
364 int i;
365
366 for (i = 0; i < length; i++)
367 value_divexact(p2[i], p1[i], lambda);
368 } /* Vector_AntiScale */
369
370 /*
371 * Puts negative of 'p1' in 'p2'
372 */
373 void Vector_Oppose(Value *p1, Value *p2, unsigned len)
374 {
375 unsigned i;
376
377 for (i = 0; i < len; ++i)
378 value_oppose(p2[i], p1[i]);
379 }
380
381 /*
382 * Return the inner product of the two Vectors 'p1' and 'p2'
383 */
384 void Inner_Product(Value *p1, Value *p2, unsigned length, Value *ip)
385 {
386 int i;
387
388 if (length != 0)
389 value_multiply(*ip, p1[0], p2[0]);
390 else
391 value_set_si(*ip, 0);
392 for(i = 1; i < length; i++)
393 value_addmul(*ip, p1[i], p2[i]);
394 } /* Inner_Product */
395
396 /*
397 * Return the maximum of the components of 'p'
398 */
399 void Vector_Max(Value *p,unsigned length, Value *max) {
400
401 Value *cp;
402 int i;
403
404 cp=p;
405 value_assign(*max,*cp);
406 cp++;
407 for (i=1;i<length;i++) {
408 value_maximum(*max,*max,*cp);
409 cp++;
410 }
411 } /* Vector_Max */
412
413 /*
414 * Return the minimum of the components of Vector 'p'
415 */
416 void Vector_Min(Value *p,unsigned length,Value *min) {
417
418 Value *cp;
419 int i;
420
421 cp=p;
422 value_assign(*min,*cp);
423 cp++;
424 for (i=1;i<length;i++) {
425 value_minimum(*min,*min,*cp);
426 cp++;
427 }
428 return;
429 } /* Vector_Min */
430
431 /*
432 * Given Vectors 'p1' and 'p2', return Vector 'p3' = lambda * p1 + mu * p2.
433 */
434 void Vector_Combine(Value *p1, Value *p2, Value *p3, Value lambda, Value mu,
435 unsigned length)
436 {
437 Value tmp;
438 int i;
439
440 value_init(tmp);
441 for (i = 0; i < length; i++) {
442 value_multiply(tmp, lambda, p1[i]);
443 value_addmul(tmp, mu, p2[i]);
444 value_assign(p3[i], tmp);
445 }
446 value_clear(tmp);
447 return;
448 } /* Vector_Combine */
449
450 /*
451 * Return 1 if 'Vec1' equals 'Vec2', otherwise return 0
452 */
453 int Vector_Equal(Value *Vec1, Value *Vec2, unsigned n)
454 {
455 int i;
456
457 for (i = 0; i < n; ++i)
458 if (value_ne(Vec1[i], Vec2[i]))
459 return 0;
460
461 return 1;
462 } /* Vector_Equal */
463
464 /*
465 * Return the component of 'p' with minimum non-zero absolute value. 'index'
466 * points to the component index that has the minimum value. If no such value
467 * and index is found, Value 1 is returned.
468 */
469 void Vector_Min_Not_Zero(Value *p,unsigned length,int *index,Value *min)
470 {
471 Value aux;
472 int i;
473
474
475 i = First_Non_Zero(p, length);
476 if (i == -1) {
477 value_set_si(*min,1);
478 return;
479 }
480 *index = i;
481 value_absolute(*min, p[i]);
482 value_init(aux);
483 for (i = i+1; i < length; i++) {
484 if (value_zero_p(p[i]))
485 continue;
486 value_absolute(aux, p[i]);
487 if (value_lt(aux,*min)) {
488 value_assign(*min,aux);
489 *index = i;
490 }
491 }
492 value_clear(aux);
493 } /* Vector_Min_Not_Zero */
494
495 /*
496 * Return the GCD of components of Vector 'p'
497 */
498 void Vector_Gcd(Value *p,unsigned length,Value *min) {
499
500 Value *q,*cq, *cp;
501 int i, Not_Zero, Index_Min=0;
502
503 q = (Value *)malloc(length*sizeof(Value));
504
505 /* Initialize all the 'Value' variables */
506 for(i=0;i<length;i++)
507 value_init(q[i]);
508
509 /* 'cp' points to vector 'p' and cq points to vector 'q' that holds the */
510 /* absolute value of elements of vector 'p'. */
511 cp=p;
512 for (cq = q,i=0;i<length;i++) {
513 value_absolute(*cq,*cp);
514 cq++;
515 cp++;
516 }
517 do {
518 Vector_Min_Not_Zero(q,length,&Index_Min,min);
519
520 /* if (*min != 1) */
521 if (value_notone_p(*min)) {
522
523 cq=q;
524 Not_Zero=0;
525 for (i=0;i<length;i++,cq++)
526 if (i!=Index_Min) {
527
528 /* Not_Zero |= (*cq %= *min) */
529 value_modulus(*cq,*cq,*min);
530 Not_Zero |= value_notzero_p(*cq);
531 }
532 }
533 else
534 break;
535 } while (Not_Zero);
536
537 /* Clear all the 'Value' variables */
538 for(i=0;i<length;i++)
539 value_clear(q[i]);
540 free(q);
541 } /* Vector_Gcd */
count_sat(unsigned int * mf)542
543 /*
544 * Given vectors 'p1', 'p2', and a pointer to a function returning 'Value type,
545 * compute p3[i] = f(p1[i],p2[i]).
546 */
547 void Vector_Map(Value *p1,Value *p2,Value *p3,unsigned length,
548 Value *(*f)(Value,Value))
549 {
550 Value *cp1, *cp2, *cp3;
551 int i;
552
553 cp1=p1;
554 cp2=p2;
555 cp3=p3;
556 for(i=0;i<length;i++) {
557 value_assign(*cp3,*(*f)(*cp1, *cp2));
558 cp1++; cp2++; cp3++;
559 }
560 return;
561 } /* Vector_Map */
562
563 /*
564 * Reduce a vector by dividing it by GCD. There is no restriction on
565 * components of Vector 'p'. Making the last element positive is *not* OK
566 * for equalities.
567 */
568 void Vector_Normalize(Value *p,unsigned length) {
569
570 Value gcd;
571 int i;
572
573 value_init(gcd);
574
575 Vector_Gcd(p,length,&gcd);
576
577 if (value_notone_p(gcd))
578 Vector_AntiScale(p, p, gcd, length);
579
580 value_clear(gcd);
581 } /* Vector_Normalize */
582
583 /*
584 * Reduce a vector by dividing it by GCD and making sure its pos-th
585 * element is positive.
586 */
587 void Vector_Normalize_Positive(Value *p,int length,int pos) {
588
589 Value gcd;
590 int i;
591
592 value_init(gcd);
593 Vector_Gcd(p,length,&gcd);
594 if (value_neg_p(p[pos]))
595 value_oppose(gcd,gcd);
596 if(value_notone_p(gcd))
597 Vector_AntiScale(p, p, gcd, length);
598 value_clear(gcd);
599 } /* Vector_Normalize_Positive */
600
601 /*
602 * Reduce 'p' by operating binary function on its components successively
603 */
604 void Vector_Reduce(Value *p,unsigned length,void(*f)(Value,Value *),Value *r) {
605
606 Value *cp;
607 int i;
608
609 cp = p;
610 value_assign(*r,*cp);
611 for(i=1;i<length;i++) {
612 cp++;
613 (*f)(*cp,r);
614 }
scan_m_face(int pos,int nb_un,Polyhedron * D,unsigned int * mf)615 } /* Vector_Reduce */
616
617 /*
618 * Sort the components of a Vector 'vector' using Heap Sort.
619 */
620 void Vector_Sort(Value *vector,unsigned n) {
621
622 int i, j;
623 Value temp;
624 Value *current_node=(Value *)0;
625 Value *left_son,*right_son;
626
627 value_init(temp);
628
629 for (i=(n-1)/2;i>=0;i--) {
630
631 /* Phase 1 : build the heap */
632 j=i;
633 value_assign(temp,*(vector+i));
634
635 /* While not a leaf */
636 while (j<=(n-1)/2) {
637 current_node = vector+j;
638 left_son = vector+(j<<1)+1;
639
640 /* If only one son */
641 if ((j<<1)+2>=n) {
642 if (value_lt(temp,*left_son)) {
643 value_assign(*current_node,*left_son);
644 j=(j<<1)+1;
645 }
646 else
647 break;
648 }
649 else {
650
651 /* If two sons */
652 right_son=left_son+1;
653 if (value_lt(*right_son,*left_son)) {
654 if (value_lt(temp,*left_son)) {
655 value_assign(*current_node,*left_son);
656 j=(j<<1)+1;
657 }
658 else
659 break;
660 }
661 else {
662 if (value_lt(temp,*right_son)) {
663 value_assign(*current_node,*right_son );
664 j=(j<<1)+2;
665 }
666 else
667 break;
668 }
669 }
670 }
671 value_assign(*current_node,temp);
672 }
673 for(i=n-1;i>0;i--) {
674
675 /* Phase 2 : sort the heap */
676 value_assign(temp, *(vector+i));
677 value_assign(*(vector+i),*vector);
678 j=0;
679
680 /* While not a leaf */
681 while (j<i/2) {
682 current_node=vector+j;
683 left_son=vector+(j<<1)+1;
684
685 /* If only one son */
686 if ((j<<1)+2>=i) {
687 if (value_lt(temp,*left_son)) {
688 value_assign(*current_node,*left_son);
689 j=(j<<1)+1;
690 }
691 else
692 break;
693 }
694 else {
695
696 /* If two sons */
697 right_son=left_son+1;
698 if (value_lt(*right_son,*left_son)) {
699 if (value_lt(temp,*left_son)) {
700 value_assign(*current_node,*left_son);
701 j=(j<<1)+1;
702 }
703 else
704 break;
705 }
706 else {
707 if (value_lt(temp,*right_son)) {
708 value_assign(*current_node,*right_son );
709 j=(j<<1)+2;
710 }
711 else
712 break;
713 }
714 }
715 }
716 value_assign(*current_node,temp);
717 }
718 value_clear(temp);
719 return;
720 } /* Vector_Sort */
721
722 /*
723 * Replaces constraint a x >= c by x >= ceil(c/a)
724 * where "a" is a common factor in the coefficients
725 * old is the constraint; v points to an initialized
726 * value that this procedure can use.
727 * Return non-zero if something changed.
728 * Result is placed in newp.
729 */
730 int ConstraintSimplify(Value *old, Value *newp, int len, Value* v)
731 {
732 /* first remove common factor of all coefficients (including "c") */
733 Vector_Gcd(old+1, len - 1, v);
734 if (value_notone_p(*v))
735 Vector_AntiScale(old+1, newp+1, *v, len-1);
736
737 Vector_Gcd(old+1, len - 2, v);
738
739 if (value_one_p(*v))
740 return 0;
741
742 Vector_AntiScale(old+1, newp+1, *v, len-2);
743 value_pdivision(newp[len-1], old[len-1], *v);
744 return 1;
745 }
746
747 int Vector_IsZero(Value * v, unsigned length) {
748 unsigned i;
749 if (value_notzero_p(v[0])) return 0;
750 else {
751 value_set_si(v[0], 1);
752 for (i=length-1; value_zero_p(v[i]); i--);
753 value_set_si(v[0], 0);
754 return (i==0);
755 }
Poly2Sat(Polyhedron * Pol,unsigned int ** L)756 }
757
758 #define MAX_CACHE_SIZE 20
759 static struct {
760 Value *p;
761 int size;
762 } cache[MAX_CACHE_SIZE];
763 static int cache_size = 0;
764
765 Value* value_alloc(int want, int *got)
766 {
767 int i;
768 Value *p;
769
770 if (cache_size) {
771 int best;
772 for (i = 0; i < cache_size; ++i) {
773 if (cache[i].size >= want) {
774 Value *p = cache[i].p;
775 *got = cache[i].size;
776 if (--cache_size != i)
777 cache[i] = cache[cache_size];
778 Vector_Set(p, 0, want);
779 return p;
780 }
781 if (i == 0)
782 best = 0;
783 else if (cache[i].size > cache[best].size)
784 best = i;
785 }
786
787 p = (Value *)realloc(cache[best].p, want * sizeof(Value));
788 *got = cache[best].size;
789 if (--cache_size != best)
790 cache[best] = cache[cache_size];
791 Vector_Set(p, 0, *got);
792 } else {
793 p = (Value *)malloc(want * sizeof(Value));
794 *got = 0;
795 }
796
797 if (!p)
798 return p;
799
800 for (i = *got; i < want; ++i)
801 value_init(p[i]);
802 *got = want;
803
804 return p;
805 }
806
807 void value_free(Value *p, int size)
808 {
GenParamPolyhedron(Polyhedron * Pol,Matrix * Rays)809 int i;
810
811 if (cache_size < MAX_CACHE_SIZE) {
812 cache[cache_size].p = p;
813 cache[cache_size].size = size;
814 ++cache_size;
815 return;
816 }
817
818 for (i=0; i < size; i++)
819 value_clear(p[i]);
820 free(p);
821 }
822
823