1 /* domains.c
2 COPYRIGHT
3 Both this software and its documentation are
4
5 Copyrighted 1997 by Vincent Loechner.
6
7 Permission is granted to copy, use, and distribute
8 for any commercial or noncommercial purpose under the terms
9 of the GNU General Public license, version 2, June 1991
10 (see file : LICENSING).
11 */
12
13 #include <stdio.h>
14 #include <math.h>
15 #include <string.h>
16 #include <stdlib.h>
17
18 #include <polylib/polylib64.h>
19
20 #include "domains.h"
21 #include "gtk_ddraw.h"
22 #include "gtk_domain.h"
23 #include "gtk_properties.h"
24 #include "gtk_context.h"
25 #include "gtk_windows.h"
26 #include "gtk_repere.h"
27 #include "visutypes.h"
28
29 int contains_int_point (Polyhedron *);
30
pgcd(int a,int b)31 int pgcd (int a, int b)
32 {
33 int r;
34
35 if (a == 0)
36 return (abs (b));
37
38 if (b == 0)
39 return (abs (a));
40
41 do
42 {
43 r = a % b;
44 a = b;
45 b = r;
46 }
47 while (r != 0);
48
49 return (abs (a));
50 }
51
52
53
54 /* Parcours de polyedre */
55 /* Usage :
56 mat = (Matrix *)Scanning_Matrix( Polyhedron *P, Matrix *order, int ws )
57 P is the Polyhedron to be scanned
58 ws the working space for the polylib calls
59 order a matrix containing the indices to be scanned (destroyed by this call) :
60 exemple :
61 order = [ 0 1 0 0 0 ]
62 [ 0 0 1 0 0 ]
63 [ 0 0 0 1 0 ]
64 to scan a polyhedron [i,j,k], with one parameter M
65
66 This call to Scanning_Matrix returns a matrix of the following form :
67 [0 i j k M cte]
68 0/1 *(!=0) 0 0 * *
69 ...
70 0/1 *!0 0 0 * *
71 0/1 * *!0 0 * *
72 ...
73 0/1 * *!0 0 * *
74 0/1 * * *!0 * *
75 ...
76 0/1 * * *!0 * *
77 */
78
scan_r(Polyhedron * P,Matrix * forall_ray,Matrix * C,int ws)79 static void scan_r (Polyhedron * P, Matrix * forall_ray, Matrix * C, int ws)
80 {
81 int l, c;
82
83 /* elimine les contraintes sur les termes restants
84 en ajoutant des lines dans ces directions. */
85 /* forall_ray = 0 1 c d d+1
86 0 0..0 1 0..0 0
87 ........
88 */
89
90 for (c = 1; c < forall_ray->NbColumns; ++c)
91 if (forall_ray->p[0][c])
92 break;
93 /* elimine la premiere ligne de la matrice forall_ray */
94 if (--forall_ray->NbRows)
95 {
96 Polyhedron *Psimp, *DU;
97 memcpy (&forall_ray->p[0][0], &forall_ray->p[1][0],
98 sizeof (Value) * forall_ray->NbColumns * forall_ray->NbRows);
99 DU = DomainAddRays (P, forall_ray, ws);
100
101 /* la borne min et max en fct des precedents est donnee par DU */
102 /* ajoute la matrice de contraintes de DU a C pour les lignes de
103 DU contenant un !=0 a l'indice c */
104
105 for (l = 0; l < DU->NbConstraints; ++l)
106 {
107 if (DU->Constraint[l][c])
108 {
109 memcpy (&C->p[C->NbRows][0], &DU->Constraint[l][0],
110 sizeof (Value) * (DU->Dimension + 2));
111 ++C->NbRows;
112 }
113 }
114
115 /* simplifie P en donnant DU comme contexte */
116 Psimp = DomainSimplify (P, DU, ws);
117 Polyhedron_Free (DU);
118
119 scan_r (Psimp, forall_ray, C, ws);
120
121 Domain_Free (Psimp);
122 }
123 else
124 {
125 /* ajoute la matrice de contraintes restant (dans P) a C */
126 for (l = 0; l < P->NbConstraints; ++l)
127 memcpy (&C->p[l + C->NbRows][0], &P->Constraint[l][0],
128 sizeof (Value) * (P->Dimension + 2));
129 C->NbRows += P->NbConstraints;
130 }
131 }
132
133 /* passer en parametre l'ordre dans lequel on veut executer la boucle. */
134 /* sous forme de vecteur ( {1,2,3,...} si ordre i->j->...) */
Scanning_Matrix(Polyhedron * P,Matrix * order,int ws)135 Matrix *Scanning_Matrix (Polyhedron * P, Matrix * order, int ws)
136 {
137 Matrix *MRes;
138
139 if (order->NbColumns != P->Dimension + 2)
140 {
141 Warn ("? Scanning_Matrix: operation on different dimensions\n");
142 return NULL;
143 }
144
145 /* if( P->next )
146 {
147 Warn("? Scanning_Matrix: can't scan a union of domains\n");
148 return NULL;
149 }
150 */
151
152
153 /* Mtmp = Matrix_Alloc( P->Dimension, P->Dimension+2 );
154 * for( l=0 ; l<P->Dimension ; ++l )
155 * {
156 * for( i=0 ; i<Mtmp->NbColumns ; ++i )
157 * Mtmp->p[l][i] = 0;
158 * if( order->p[l]<1 || order->p[l]>P->Dimension )
159 * fprintf(stderr,"? Scanning_Matrix: incorrect order vector\n");
160 * Mtmp->p[l][order->p[l]] = 1;
161 * }
162 */
163
164 MRes = Matrix_Alloc (ws, P->Dimension + 2);
165 MRes->NbRows = 0;
166
167 scan_r (P, order, MRes, ws);
168
169 /* Matrix_Free( Mtmp ); */
170
171 if (MRes->NbRows == 0)
172 {
173 Matrix_Free (MRes);
174 return NULL;
175 }
176 else
177 { /* trop de lignes sont allouees...
178 on recopie cette matrice dans une nouvelle */
179 Matrix *M;
180 M = Matrix_Alloc (MRes->NbRows, MRes->NbColumns);
181 memcpy (&M->p[0][0], &MRes->p[0][0],
182 sizeof (Value) * MRes->NbRows * MRes->NbColumns);
183 Matrix_Free (MRes);
184 return (M);
185 }
186 }
187
188
189
190 /* parcours d'un polyedre convexe (pas union) avec appel de fonction */
191 /* pour chaque pt entier */
192 /* USAGE : FORALL_P_do( Polyhedron *P, int (*f)(Value *) )
193 la fonction f est appelee en passant le pt courant en parametre */
194 /* le parcours s'arrete si f renvoie FALSE */
195 /* la fonction renvoie FALSE si elle s'est arretee avant la fin */
196 /* c = indice parcourant 0 -> NbC-2 inclus */
197 static int
FORALL_P_at_pt_do(Matrix * C,int c,Value * pt,int (* f)(Value *))198 FORALL_P_at_pt_do (Matrix * C, int c, Value * pt, int (*f) (Value *))
199 {
200 int l, i, j;
201 int min_num = 0, max_num = 0, new_num;
202 int min_den = 1, max_den = 1, new_den;
203 int trouve_min, trouve_max;
204 int g, den;
205
206 /*
207 Cherche une ligne du type :
208 0 1 ... c+1 c+2 .. NbC-2 NbC-1
209 C = * * .... * 0 .... 0 *
210 */
211
212 trouve_min = trouve_max = 0;
213
214 /* on parcourt tous les !=0 sur cette colonne */
215 for (l = 0; l < C->NbRows; ++l)
216 {
217 if (C->p[l][c + 1] == 0)
218 continue;
219
220 /* si sur cette ligne,
221 il y a un !=0 apres cet indice on ignore la ligne */
222 for (j = c + 2; j < C->NbColumns - 1; ++j)
223 if (C->p[l][j])
224 break;
225 if (j < C->NbColumns - 1)
226 continue;
227
228 /* OK, la ligne 'l' contient bien une ref utile a cette variable */
229
230 /* cherche la valeur de i : cte +val des ind. precedent */
231 new_num = C->p[l][C->NbColumns - 1];
232 for (j = 1; j <= c; ++j)
233 new_num += pt[j - 1] * C->p[l][j];
234 new_den = C->p[l][c + 1];
235 if (new_den < 0) /* denominateurs toujours > 0 */
236 new_den = -new_den;
237 else
238 new_num = -new_num;
239 if (C->p[l][0])
240 {
241 /* inegalite */
242 if (C->p[l][c + 1] > 0)
243 {
244 /* new est un nouveau min */
245 if (trouve_min)
246 {
247 if (new_num * min_den > min_num * new_den)
248 {
249 min_num = new_num;
250 min_den = new_den;
251 }
252 }
253 else
254 {
255 min_num = new_num;
256 min_den = new_den;
257 trouve_min = 1;
258 }
259 }
260 else
261 { /* new est un nouveau max */
262 if (trouve_max)
263 {
264 if (new_num * max_den < max_num * new_den)
265 {
266 max_num = new_num;
267 max_den = new_den;
268 }
269 }
270 else
271 {
272 max_num = new_num;
273 max_den = new_den;
274 trouve_max = 1;
275 }
276 }
277 /* si le min est devenu > au max pas la peine de continuer */
278 if (trouve_min && trouve_max)
279 if (max_num * min_den < min_num * max_den)
280 return 1;
281 }
282 else
283 { /* egalite */
284 /* si on a deja un min ou un max qui empiete on s'arrete la. */
285 if (trouve_min && (new_num * min_den < min_num * new_den))
286 return 1;
287 if (trouve_max && (new_num * max_den > max_num * new_den))
288 return 1;
289 trouve_min = 1;
290 trouve_max = 1;
291 min_num = max_num = new_num;
292 min_den = max_den = new_den;
293 }
294 }
295
296 if (!trouve_min || !trouve_max)
297 {
298 fprintf (stderr, "Error in scan polyhedron : non closed polyhedron\n");
299 Matrix_Print (stderr, "%3d ", C);
300
301 Warn ("Can't visualize an infinite polyhedron\n");
302
303 return 0;
304 }
305
306
307 g = pgcd (min_den, max_den);
308 den = min_den * max_den / g;
309 min_num = min_num * max_den / g;
310 max_num = max_num * min_den / g;
311
312 /* se place sur le premier "vrai" */
313 if (min_num < 0)
314 min_num += (-min_num) % den;
315 else if ((min_num) % den)
316 min_num += den - (min_num) % den;
317
318 /* decrit la boucle */
319 if (c == C->NbColumns - 3)
320 {
321 for (i = min_num; i <= max_num; i += den)
322 {
323 pt[c] = i / den;
324 if (!(*f) (pt))
325 return 0;
326 }
327 }
328 else
329 {
330 for (i = min_num; i <= max_num; i += den)
331 {
332 pt[c] = i / den;
333 if (!FORALL_P_at_pt_do (C, c + 1, pt, f))
334 return 0;
335 }
336 }
337 return 1;
338 }
339
FORALL_P_do(Polyhedron * P,Value * pt,int (* f)(Value *))340 int FORALL_P_do (Polyhedron * P, Value * pt, int (*f) (Value *))
341 {
342 Matrix *constraints;
343 Matrix *o;
344 int i;
345
346 if (!P || P->NbRays == 0)
347 return 1;
348
349 o = Matrix_Alloc (P->Dimension, P->Dimension + 2);
350 memset (&o->p[0][0], 0, o->NbRows * o->NbColumns * sizeof (Value));
351 for (i = 0; i < o->NbRows; ++i)
352 o->p[i][i + 1] = 1;
353
354 constraints = Scanning_Matrix (P, o, WS);
355 Matrix_Free (o);
356
357 if (constraints)
358 return (FORALL_P_at_pt_do (constraints, 0, pt, f));
359 else
360 return 1;
361
362
363 }
364
365
_CIP_ok(Value * pt)366 static int _CIP_ok (Value * pt)
367 {
368 return 0;
369 }
370
contains_int_point(Polyhedron * P)371 int contains_int_point (Polyhedron * P)
372 {
373 int r;
374 Value *pt;
375 pt = malloc (sizeof (Value) * P->Dimension);
376 r = FORALL_P_do (P, pt, &_CIP_ok);
377 free (pt);
378 return (!r);
379 }
380
381 static int compt;
compte_nb_pts_poly(Value * pt)382 static int compte_nb_pts_poly (Value * pt)
383 {
384 ++compt;
385 return 1;
386 }
387
388 /* compte le nb de pts entiers dans un polyedre */
nb_pts_poly(Polyhedron * P)389 int nb_pts_poly (Polyhedron * P)
390 {
391 Value *pt;
392 pt = malloc (sizeof (Value) * P->Dimension);
393 compt = 0;
394 FORALL_P_do (P, pt, &compte_nb_pts_poly);
395 free (pt);
396
397 return (compt);
398 }
399
400 /* teste si le point de coord. pt est dans P */
DansPolyedre(Polyhedron * P,Value * pt)401 int DansPolyedre (Polyhedron * P, Value * pt)
402 {
403 int c, i, r;
404
405 for (c = 0; c < P->NbConstraints; ++c)
406 {
407 r = 0;
408 for (i = 0; i < P->Dimension; ++i)
409 {
410 r += pt[i] * P->Constraint[c][i + 1];
411 }
412 r += P->Constraint[c][P->Dimension + 1]; /* constante */
413
414 if (P->Constraint[c][0]) /* inegalite ou egalite */
415 {
416 if (r < 0)
417 return (0);
418 }
419 else
420 {
421 if (r != 0)
422 return (0);
423 }
424 }
425 return (1);
426 }
427
428 /* teste si le point de coord. pt est dans le domaine P (convexe ou non) */
DansDomaine(Polyhedron * P,Value * pt)429 gboolean DansDomaine (Polyhedron * P, Value * pt)
430 {
431 Polyhedron * tmp;
432
433 for (tmp = P; tmp; tmp = tmp->next)
434 if( DansPolyedre (tmp, pt) )
435 return( TRUE );
436 return( FALSE );
437 }
438
439
440 /************************************************************************/
441 /* NP_Domain() */
442 /************************************************************************/
443 /* calcul du polyhedre non parametre correspondant a l'instantiation
444 des params dans p. */
NP_Domain(Polyhedron * P)445 static Polyhedron *NP_Domain (Polyhedron * P)
446 {
447 Matrix *CNP; /* Non-Parameterized constraints */
448 Polyhedron *DNP;
449 int l, k;
450
451 if (!P)
452 return (NULL);
453
454 CNP = Matrix_Alloc (P->NbConstraints, DimNP + 2);
455 for (l = 0; l < P->NbConstraints; ++l)
456 {
457 /* les dnp+1 premieres colonnes */
458 for (k = 0; k <= DimNP; ++k)
459 CNP->p[l][k] = P->Constraint[l][k];
460
461 /* la constante : */
462 CNP->p[l][DimNP + 1] = P->Constraint[l][1 + P->Dimension];
463 for (k = 0; k < NbParam; ++k)
464 CNP->p[l][DimNP + 1] +=
465 ValeurParam[k] * P->Constraint[l][DimNP + 1 + k];
466 }
467 DNP = Constraints2Polyhedron (CNP, WS);
468 Matrix_Free (CNP);
469
470 DNP = AddPolyToDomain( DNP, NP_Domain(P->next) );
471 if( DNP )
472 if( emptyQ(DNP) )
473 {
474 Domain_Free( DNP );
475 return( NULL );
476 }
477 return( DNP );
478
479
480 }
481
482
483
clear_A_domains()484 void clear_A_domains ()
485 {
486 if (U_NPDomains)
487 Domain_Free (U_NPDomains);
488 if (A_Union)
489 Domain_Free (A_Union);
490 if (A_Conv)
491 Domain_Free (A_Conv);
492 if (U_IntPDomains)
493 Domain_Free (U_IntPDomains);
494
495 A_Union = NULL;
496 A_Conv = NULL;
497 U_NPDomains = NULL;
498 U_IntPDomains = NULL;
499 }
500
501 /* Cr�e l'union des domaines que l'on veut afficher pour pouvoir calculer
502 les extr�mas de la figure � visualiser... */
compute_A_domains()503 void compute_A_domains ()
504 {
505 Domain *d;
506
507 if (U_NPDomains)
508 Domain_Free (U_NPDomains);
509 /* Construction de l'union des domaines (instanci�s) � afficher */
510 U_NPDomains = NULL;
511 for( d = DomainList ; d ; d = d->next )
512 {
513 if (d->visual != NULL)
514 {
515 /* U_NPDomains = DomainUnion (d->visual->NPdomain, U_NPDomains, WS); */
516 Polyhedron *p;
517 for( p=d->visual->NPdomain ; p ; p=p->next )
518 U_NPDomains = AddPolyToDomain(Polyhedron_Copy(p), U_NPDomains);
519 }
520 }
521
522
523 /* Cr�e l'enveloppe convexe de l'union des domaines que l'on veut afficher pour
524 calculer les extr�mas de la figure � visualiser... */
525 if (A_Conv)
526 Domain_Free (A_Conv);
527
528 A_Conv = DomainConvex( U_NPDomains, WS);
529
530 /********************** limite le bord � -10,10 dans toutes les dimensions **************/
531 /*
532 r->NbRows = A_Conv->Dimension*2;
533 for( j=0 ; j<A_Conv->Dimension ; j++ )
534 {
535 for( i=1 ; i<DimNP+1 ; i++ )
536 {
537 r->p[2*j][i] = 0;
538 r->p[2*j+1][i] = 0;
539 }
540 r->p[2*j][0] = 1;
541 r->p[2*j][j+1] = -1;
542 r->p[2*j][DimNP + 1] = 10;
543
544 r->p[2*j+1][0] = 1;
545 r->p[2*j+1][j+1] = 1;
546 r->p[2*j+1][DimNP + 1] = 10;
547 }
548 A_Conv = DomainAddConstraints( A_Conv, r, WS );
549
550 Matrix_Free (r);
551 */
552
553 return;
554 }
555
556 /* Cr�e l'union des domaines dont les points entiers sont visualis�s... */
union_intpoints_domains()557 void union_intpoints_domains ()
558 {
559 Domain *d;
560
561 if (U_IntPDomains)
562 Domain_Free (U_IntPDomains);
563
564 U_IntPDomains = NULL;
565 for( d = DomainList ; d ; d = d->next )
566 {
567 if (d->visual && d->visual->intpoints )
568 {
569 /* U_IntPDomains = DomainUnion (d->visual->NPdomain, U_IntPDomains, WS); */
570 Polyhedron *p;
571 for( p=d->visual->NPdomain ; p ; p=p->next )
572 U_IntPDomains = AddPolyToDomain(Polyhedron_Copy(p), U_IntPDomains);
573 }
574 }
575
576 return;
577 }
578
579 /* Proc�dure d'affichage de domaines param�tr�s */
Affich_P_domains(Domain * domain,gint row)580 void Affich_P_domains (Domain * domain, gint row)
581 {
582 Polyhedron *tmp_domain;
583 gint Rrow;
584 gboolean delete;
585 char *Status;
586
587 delete = FALSE;
588 if (ContextOK == TRUE)
589 {
590 /* Les param�tres et le contexte sont OK */
591
592 if (domain->Pdomain->Dimension - NbParam != DimNP)
593 {
594 /* La dimension du domaine n'est pas coh�rente avec les
595 initialisations de DimNP et du contexte (Parano�a: �a
596 ne peut pas arriver car dans ce cas le domaine n'est
597 pas dans la liste graphique!) */
598 printf
599 ("Domain %s can't be vizualized with this context!\n",
600 domain->name);
601 delete = TRUE;
602 }
603 else
604 {
605 /* Le domaine aura la bonne dimension apr�s instanciation */
606 tmp_domain = NP_Domain (domain->Pdomain);
607 if( tmp_domain )
608 {
609 // A CHANGER 0 &
610 if ( value_zero_p
611 (tmp_domain->Ray[0][tmp_domain->Dimension + 1]))
612 {
613 /* Le domaine � afficher n'est pas born� dans ce contexte */
614 printf("Domain %s is not bounded!\n", domain->name);
615 delete = TRUE;
616 }
617 else
618 {
619 /* Le domaine instanci� est born� */
620 /* Modifier la structure du domaine en cons�quence */
621 if (!domain->visual)
622 {
623 /* Ce domaine n'�tait pas affich� */
624 domain->visual =
625 (VisualDomain *) malloc (sizeof (VisualDomain));
626 // A CHANGER domain->visual->intpoints = FALSE;
627 domain->visual->intpoints = TRUE;
628 domain->visual->Lpoints = NULL;
629 domain->visual->NPdomain = NULL;
630 /* M-�-j de la liste graphique */
631 gtk_clist_set_foreground (GTK_CLIST (clist), row,
632 domain->color);
633 gtk_clist_set_text (GTK_CLIST (clist), row, 1, "Visu");
634 }
635 if (domain->visual->NPdomain)
636 Domain_Free (domain->visual->NPdomain);
637 domain->visual->NPdomain = tmp_domain;
638 if (domain->visual->Lpoints)
639 ListValuePoints_Free (domain->visual->Lpoints);
640 domain->visual->Lpoints = NULL;
641 domain->visual->Lpoints = ListValuePoints_Update (domain);
642 }
643 }
644 else
645 {
646 printf ("Empty domain %s!\n", domain->name);
647 delete = TRUE;
648 }
649
650 }
651 }
652 else
653 {
654 /* ContextOK = FALSE */
655 /* pas grave !
656 Warn ("Bad parameters!");
657 delete = TRUE;
658 */
659 }
660
661 if (delete)
662 {
663 if (domain->visual)
664 {
665 /* Ce domaine �tait affich� avant le changement de param�tres, il
666 faut le masquer... */
667 Rrow = gtk_clist_find_row_from_data (GTK_CLIST (clist), domain);
668 if (Rrow == -1)
669 {
670 Warn ("Internal problem: to mask a domain!");
671 return;
672 }
673 else
674 {
675 delete = gtk_clist_get_text (GTK_CLIST (clist), Rrow, 1, &Status);
676 if (strcmp (Status, "Visu") == 0)
677 gtk_clist_select_row (GTK_CLIST (clist), Rrow, 1);
678 }
679 }
680 }
681 compute_A_domains ();
682
683 calc_rep_extr ();
684 return;
685 }
686
Affich_NOT_P_domains(Domain * domain,gint row)687 void Affich_NOT_P_domains (Domain * domain, gint row)
688 {
689 if (domain->Pdomain->Dimension != DimNP)
690 /* La dimension du domaine n'est pas coh�rente avec les
691 initialisations de DimNP et du contexte */
692 printf
693 ("Domain %s can't be visualized!\n", domain->name);
694 else
695 {
696 /* Le domaine a la bonne dimension */
697 if( domain->Pdomain && domain->Pdomain->NbRays )
698 {
699 if (value_zero_p
700 (domain->Pdomain->Ray[0][domain->Pdomain->Dimension + 1]))
701 /* Le domaine � afficher n'est pas born� dans ce contexte */
702 printf
703 ("Domain %s is not bounded!\n", domain->name);
704 else
705 {
706 /* Le domaine est born� */
707
708 /* Modifier la structure du domaine en cons�quence */
709 domain->visual = (VisualDomain *) malloc (sizeof (VisualDomain));
710 domain->visual->intpoints = TRUE;
711 domain->visual->Lpoints = NULL;
712 domain->visual->NPdomain = Domain_Copy (domain->Pdomain);
713 /* M-�-j de la liste graphique */
714 gtk_clist_set_foreground (GTK_CLIST (clist), row, domain->color);
715 gtk_clist_set_text (GTK_CLIST (clist), row, 1, "Visu");
716 }
717 }
718 else
719 printf ("Empty domain %s!\n", domain->name);
720
721 }
722 /* Il y a des domaines � afficher */
723 compute_A_domains ();
724
725 calc_rep_extr ();
726 return;
727 }
728