source/piplib.c

   1 /******************************************************************************
   2  *                     PIP : Parametric Integer Programming                   *
   3  ******************************************************************************
   4  *                                 piplib.c                                   *
   5  ******************************************************************************
   6  *                                                                            *
   7  * Copyright Paul Feautrier, 1988-2005                                        *
   8  *                                                                            *
   9  * This is free software; you can redistribute it and/or modify it under the  *
  10  * terms of the GNU General Public License as published by the Free Software  *
  11  * Foundation; either version 2 of the License, or (at your option) any later *
  12  * version.                                                                   *
  13  *                                                                            *
  14  * This software is distributed in the hope that it will be useful, but       *
  15  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
  16  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
  17  * for more details.                                                          *
  18  *                                                                            *
  19  * You should have received a copy of the GNU General Public License along    *
  20  * with software; if not, write to the Free Software Foundation, Inc.,        *
  21  * 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA                     *
  22  *                                                                            *
  23  * Written by Cedric Bastoul                                                  *
  24  *                                                                            *
  25  ******************************************************************************/
  26
  27 /* Premiere version du 30 juillet 2001. */
  28
  29 # include <stdlib.h>
  30 # include <stdio.h>
  31 # include <ctype.h>
  32 #include <string.h>
  33
  34 #include "pip.h"
  35 #define min(x,y) ((x) < (y)? (x) : (y))
  36
  37 Entier UN;
  38 Entier ZERO;
  39
  40 long int cross_product, limit;
  41 int allocation, comptage;
  42 int verbose = 0;
  43 int profondeur = 0;
  44 int compa_count;
  45 int deepest_cut = 0;
  46
  47 FILE *dump = NULL;
  48 char dump_name[] = "PipXXXXXX";
  49
  50 /* Larger line buffer to accomodate Frédo Vivien exemples. A version
  51 handling arbitrary line length should be written ASAP.
  52 */
  53
  54 #define INLENGTH 2048
  55
  56 char inbuff[INLENGTH];
  57 int inptr = 256;
  58 int proviso = 0;
  59
  60
  61 /******************************************************************************
  62  *                 Fonctions d'acquisition de données (ex-maind.c)            *
  63  ******************************************************************************/
  64
  65
  66 int dgetc(FILE *foo)
  67 {
  68  char *p;
  69  if(inptr >= proviso)
  70    {p = fgets(inbuff, INLENGTH, foo);
  71     if(p == NULL) return EOF;
  72     proviso = min(INLENGTH, strlen(inbuff));
  73     inptr = 0;
  74     if(verbose > 2) fprintf(dump, "-- %s", inbuff);
  75   }
  76  return inbuff[inptr++];
  77 }
  78
  79
  80 #if defined(LINEAR_VALUE_IS_MP)
  81 int dscanf(FILE *foo, Entier  val)
  82 #else
  83 int dscanf(FILE *foo, Entier *val)
  84 #endif
  85 {
  86  char * p;
  87  int c;
  88
  89  for(;inptr < proviso; inptr++)
  90    if(inbuff[inptr] != ' ' && inbuff[inptr] != '\n' && inbuff[inptr] != '\t')
  91                                 break;
  92  while(inptr >= proviso)
  93    {p = fgets(inbuff, 256, foo);
  94     if(p == NULL) return EOF;
  95     proviso = strlen(inbuff);
  96     if(verbose > 2) {
  97       fprintf(dump, ".. %s", inbuff);
  98       fflush(dump);
  99     }
 100     for(inptr = 0; inptr < proviso; inptr++)
 101        if(inbuff[inptr] != ' '
 102        && inbuff[inptr] != '\n'
 103        && inbuff[inptr] != '\t') break;
 104   }
 105  #if defined(LINEAR_VALUE_IS_MP)
 106  if(gmp_sscanf(inbuff+inptr, GMP_INPUT_FORMAT, val) != 1)
 107  #else
 108  if(sscanf(inbuff+inptr, FORMAT, val) != 1)
 109  #endif
 110  return -1;
 111
 112  for(; inptr < proviso; inptr++)
 113         if((c = inbuff[inptr]) != '-' && !isdigit(c)) break;
 114  return 0;
 115 }
 116
 117
 118 /******************************************************************************
 119  *                    Fonctions d'affichage des structures                    *
 120  ******************************************************************************/
 121
 122
 123 /* Fonction pip_matrix_print :
 124  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
 125  * que contient la structure de type PipMatrix qu'elle recoit en parametre.
 126  * Premiere version : Ced. 29 juillet 2001.
 127  */
 128 void pip_matrix_print(FILE * foo, PipMatrix * Mat)
 129 { Entier * p;
 130   int i, j ;
 131   unsigned NbRows, NbColumns ;
 132
 133   fprintf(foo,"%d %d\n", NbRows=Mat->NbRows, NbColumns=Mat->NbColumns) ;
 134   for (i=0;i<NbRows;i++)
 135   { p=*(Mat->p+i) ;
 136     for (j=0;j<NbColumns;j++)
 137     #if defined(LINEAR_VALUE_IS_MP)
 138     { fprintf(foo," ") ;
 139       mpz_out_str(foo,10,*p++) ;
 140     }
 141     #else
 142     fprintf(foo," "FORMAT, *p++) ;
 143     #endif
 144     fprintf(foo, "\n") ;
 145   }
 146 }
 147
 148
 149 /* Fonction pip_vector_print :
 150  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
 151  * que contient la structure de type PipVector qu'elle recoit en parametre.
 152  * Premiere version : Ced. 20 juillet 2001.
 153  */
 154 void pip_vector_print(FILE * foo, PipVector * vector)
 155 { int i ;
 156
 157   if (vector != NULL)
 158   { fprintf(foo,"#[") ;
 159     for (i=0;i<vector->nb_elements;i++)
 160     { fprintf(foo," ") ;
 161       #if defined(LINEAR_VALUE_IS_MP)
 162       mpz_out_str(foo,10,vector->the_vector[i]) ;
 163       if (mpz_cmp(vector->the_deno[i],UN) != 0)
 164       #else
 165       fprintf(foo,FORMAT,vector->the_vector[i]) ;
 166       if (vector->the_deno[i] != UN)
 167       #endif
 168       { fprintf(foo,"/") ;
 169         #if defined(LINEAR_VALUE_IS_MP)
 170         mpz_out_str(foo,10,vector->the_deno[i]) ;
 171         #else
 172         fprintf(foo,FORMAT,vector->the_deno[i]) ;
 173         #endif
 174       }
 175     }
 176     fprintf(foo,"]") ;
 177   }
 178 }
 179
 180
 181 /* Fonction pip_newparm_print :
 182  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
 183  * que contient la structure de type PipNewparm qu'elle recoit en parametre.
 184  * Le parametre indent est le nombre d'espaces blancs en debut de chaque
 185  * ligne avant indentation. Une valeur negative de indent signifie qu'on ne
 186  * desire pas d'indentation.
 187  * Premiere version : Ced. 18 octobre 2001.
 188  */
 189 void pip_newparm_print(FILE * foo, PipNewparm * newparm, int indent)
 190 { int i ;
 191
 192   if (newparm != NULL)
 193   { do
 194     { for (i=0;i<indent;i++) fprintf(foo," ") ;             /* Indent. */
 195       fprintf(foo,"(newparm ") ;
 196       fprintf(foo,"%d",newparm->rank) ;
 197       fprintf(foo," (div ") ;
 198       pip_vector_print(foo,newparm->vector) ;
 199       fprintf(foo," ") ;
 200       #if defined(LINEAR_VALUE_IS_MP)
 201       mpz_out_str(foo,10,newparm->deno) ;
 202       #else
 203       fprintf(foo,FORMAT,newparm->deno) ;
 204       #endif
 205       fprintf(foo,"))\n") ;
 206     }
 207     while ((newparm = newparm->next) != NULL) ;
 208   }
 209 }
 210
 211
 212 /* Fonction pip_list_print :
 213  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
 214  * que contient la structure de type PipList qu'elle recoit en parametre.
 215  * Le parametre indent est le nombre d'espaces blancs en debut de chaque
 216  * ligne avant indentation. Une valeur negative de indent signifie qu'on ne
 217  * desire pas d'indentation.
 218  * Premiere version : Ced. 18 octobre 2001.
 219  * 16 novembre 2005 : Ced. Prise en compte du cas list->vector == NULL,
 220  *                         jusque là impossible.
 221  */
 222 void pip_list_print(FILE * foo, PipList * list, int indent)
 223 { int i ;
 224
 225   if (list == NULL)
 226   { for (i=0;i<indent;i++) fprintf(foo," ") ;               /* Indent. */
 227     fprintf(foo,"()\n") ;
 228   }
 229   else
 230   { for (i=0;i<indent;i++) fprintf(foo," ") ;               /* Indent. */
 231     fprintf(foo,"(list\n") ;
 232     do
 233     { if (list->vector != NULL)
 234       { for (i=0;i<indent+1;i++) fprintf(foo," ") ;         /* Indent. */
 235         pip_vector_print(foo,list->vector) ;
 236         fprintf(foo,"\n") ;
 237       }
 238     }
 239     while ((list = list->next) != NULL) ;
 240     for (i=0;i<indent;i++) fprintf(foo," ") ;               /* Indent. */
 241     fprintf(foo,")\n") ;
 242   }
 243 }
 244
 245
 246 /* Fonction pip_quast_print :
 247  * Cette fonction se charge d'imprimer sur le flux 'foo' les informations
 248  * que contient la structure de type PipQuast qu'elle recoit en parametre.
 249  * Le parametre indent est le nombre d'espaces blancs en debut de chaque
 250  * ligne avant indentation. Une valeur negative de indent signifie qu'on ne
 251  * desire pas d'indentation.
 252  * 20 juillet 2001 : Premiere version, Ced.
 253  * 18 octobre 2001 : eclatement.
 254  */
 255 void pip_quast_print(FILE * foo, PipQuast * solution, int indent)
 256 { int i ;
 257   PipVector * vector ;
 258   int new_indent = indent >= 0 ? indent+1 : indent;
 259
 260   if (solution != NULL)
 261   { pip_newparm_print(foo,solution->newparm,indent) ;
 262     if (solution->condition == NULL) {
 263         pip_list_print(foo, solution->list, indent);
 264         /* Possible dual solution */
 265         if (solution->next_then)
 266             pip_quast_print(foo, solution->next_then, new_indent);
 267     } else
 268     { for (i=0;i<indent;i++) fprintf(foo," ") ;             /* Indent. */
 269       fprintf(foo,"(if ") ;
 270       pip_vector_print(foo,solution->condition) ;
 271       fprintf(foo,"\n") ;
 272       pip_quast_print(foo, solution->next_then, new_indent);
 273       pip_quast_print(foo, solution->next_else, new_indent);
 274       for (i=0;i<indent;i++) fprintf(foo," ") ;             /* Indent. */
 275       fprintf(foo,")\n") ;
 276     }
 277   }
 278   else
 279   { for (i=0;i<indent;i++) fprintf(foo," ") ;               /* Indent. */
 280     fprintf(foo,"void\n") ;
 281   }
 282 }
 283
 284
 285 /* Function pip_options_print:
 286  * This function prints the content of a PipOptions structure (options)
 287  * into a file (foo, possibly stdout).
 288  * March 17th 2003: first version.
 289  */
 290 void * pip_options_print(FILE * foo, PipOptions * options)
 291 { fprintf(foo,"Option setting is:\n") ;
 292   fprintf(foo,"Nq          =%d\n",options->Nq) ;
 293   fprintf(foo,"Verbose     =%d\n",options->Verbose) ;
 294   fprintf(foo,"Simplify    =%d\n",options->Simplify) ;
 295   fprintf(foo,"Deepest_cut =%d\n",options->Deepest_cut) ;
 296   fprintf(foo,"Maximize    =%d\n",options->Maximize) ;
 297   fprintf(foo,"Urs_parms   =%d\n",options->Urs_parms);
 298   fprintf(foo,"Urs_unknowns=%d\n",options->Urs_unknowns);
 299   fprintf(foo,"\n") ;
 300 }
 301
 302
 303 /******************************************************************************
 304  *                       Fonctions de liberation memoire                      *
 305  ******************************************************************************/
 306
 307
 308 /* Fonction pip_matrix_free :
 309  * Cette fonction libere la memoire reservee a la structure de type PipMatrix
 310  * que pointe son parametre.
 311  * Premiere version : Ced. 29 juillet 2001.
 312  */
 313 void pip_matrix_free(PipMatrix * matrix)
 314 {
 315   #if defined(LINEAR_VALUE_IS_MP)
 316   int i, j ;
 317   Entier * p ;
 318
 319   p = matrix->p_Init ;
 320   for (i=0;i<matrix->p_Init_size;i++)
 321   mpz_clear(*p++) ;
 322   #endif
 323
 324   if (matrix != NULL)
 325   { free(matrix->p_Init) ;
 326     free(matrix->p) ;
 327     free(matrix) ;
 328   }
 329 }
 330
 331
 332 /* Fonction pip_vector_free :
 333  * Cette fonction libere la memoire reservee a la structure de type PipVector
 334  * que pointe son parametre.
 335  * 20 juillet 2001 : Premiere version, Ced.
 336  * 18 octobre 2001 : simplification suite a l'eclatement de PipVector.
 337  * 16 novembre 2005 : Ced. Prise en compte du cas NULL.
 338  */
 339 void pip_vector_free(PipVector * vector)
 340 { int i ;
 341
 342   if (vector != NULL)
 343   {
 344     #if defined(LINEAR_VALUE_IS_MP)
 345     for (i=0;i<vector->nb_elements;i++)
 346     { mpz_clear(vector->the_vector[i]);
 347       mpz_clear(vector->the_deno[i]);
 348     }
 349     #endif
 350
 351     free(vector->the_vector) ;
 352     free(vector->the_deno) ;
 353     free(vector) ;
 354   }
 355 }
 356
 357
 358 /* Fonction pip_newparm_free :
 359  * Cette fonction libere la memoire reservee a la structure de type PipNewparm
 360  * que pointe son parametre. Sont liberes aussi tous les elements de la
 361  * liste chainee dont il pouvait etre le depart.
 362  * Premiere version : Ced. 18 octobre 2001.
 363  */
 364 void pip_newparm_free(PipNewparm * newparm)
 365 { PipNewparm * next ;
 366
 367   while (newparm != NULL)
 368   { next = newparm->next ;
 369     #if defined(LINEAR_VALUE_IS_MP)
 370     mpz_clear(newparm->deno);
 371     #endif
 372     pip_vector_free(newparm->vector) ;
 373     free(newparm) ;
 374     newparm = next ;
 375   }
 376 }
 377
 378
 379 /* Fonction pip_list_free :
 380  * Cette fonction libere la memoire reservee a la structure de type PipList
 381  * que pointe son parametre. Sont liberes aussi tous les elements de la
 382  * liste chainee dont il pouvait etre le depart.
 383  * Premiere version : Ced. 18 octobre 2001.
 384  */
 385 void pip_list_free(PipList * list)
 386 { PipList * next ;
 387
 388   while (list != NULL)
 389   { next = list->next ;
 390     pip_vector_free(list->vector) ;
 391     free(list) ;
 392     list = next ;
 393   }
 394 }
 395
 396
 397 /* Fonction pip_quast_free :
 398  * Cette fonction libere la memoire reservee a la structure de type
 399  * PipSolution que pointe son parametre. Sont liberees aussi toutes les
 400  * differentes listes chainees qui pouvaient en partir.
 401  * 20 juillet 2001 : Premiere version, Ced.
 402  * 18 octobre 2001 : simplification suite a l'eclatement de PipVector.
 403  */
 404 void pip_quast_free(PipQuast * solution)
 405 {
 406     if (!solution)
 407         return;
 408     pip_newparm_free(solution->newparm) ;
 409
 410     pip_list_free(solution->list) ;
 411
 412     pip_vector_free(solution->condition);
 413     pip_quast_free(solution->next_then);
 414     pip_quast_free(solution->next_else);
 415     free(solution) ;
 416 }
 417
 418
 419 /* Funtion pip_options_free:
 420  * This function frees the allocated memory for a PipOptions structure.
 421  * March 15th 2003: first version.
 422  */
 423 void pip_options_free(PipOptions * options)
 424 { free(options) ;
 425 }
 426
 427
 428 /******************************************************************************
 429  *                     Fonction d'initialisation des options                  *
 430  ******************************************************************************/
 431
 432
 433 /* Funtion pip_options_init:
 434  * This function allocates the memory for a PipOptions structure and fill the
 435  * options with the default values.
 436  ********
 437  * Nq est un booleen renseignant si on cherche
 438  * une solution entiere (vrai=1) ou non (faux=0). Verbose est un booleen
 439  * permettant de rendre Pip bavard (Verbose a vrai=1), il imprimera
 440  * alors la plupart de ses traitements dans le fichier dont le nom est
 441  * dans la variable d'environnement DEBUG, ou si DEBUG
 442  * n'est pas placee, dans un nouveau fichier de nom genere par mkstemp, si
 443  * Verbose est a faux=0, Pip restera muet. Simplify est un booleen permettant
 444  * de demander a Pip de simplifier sa solution (en eliminant les formes de type
 445  * 'if #[...] () ()') quand il est a vrai=1, ou non quand il est a faux=0. Max
 446  * n'est pas encore utilise et doit etre mis a 0.
 447  ********
 448  * March 15th 2003: first version.
 449  */
 450 PipOptions * pip_options_init(void)
 451 { PipOptions * options ;
 452
 453   options = (PipOptions *)malloc(sizeof(PipOptions)) ;
 454   /* Default values of the options. */
 455   options->Nq          = 1 ;  /* Integer solution. */
 456   options->Verbose     = 0 ;  /* No comments. */
 457   options->Simplify    = 0 ;  /* Do not simplify solutions. */
 458   options->Deepest_cut = 0 ;  /* Do not use deepest cut algorithm. */
 459   options->Maximize    = 0 ;  /* Do not compute maximum. */
 460   options->Urs_parms   = 0 ;  /* All parameters are non-negative. */
 461   options->Urs_unknowns= 0 ;  /* All unknows are non-negative. */
 462   options->Compute_dual= 0;   /* Don't compute dual variables. */
 463
 464   return options ;
 465 }
 466
 467
 468 /******************************************************************************
 469  *                     Fonctions d'acquisition de matrices                    *
 470  ******************************************************************************/
 471
 472
 473 /* Fonction pip_matrix_alloc :
 474  * Fonction (tres) modifiee de Matrix_Alloc de la polylib. Elle alloue l'espace
 475  * memoire necessaire pour recevoir le contenu d'une matrice de NbRows lignes
 476  * et de NbColumns colonnes, et initialise les valeurs a 0. Elle retourne un
 477  * pointeur sur l'espace memoire alloue.
 478  * Premiere version : Ced. 18 octobre 2001.
 479  */
 480 PipMatrix * pip_matrix_alloc(unsigned NbRows, unsigned NbColumns)
 481 { PipMatrix * matrix ;
 482   Entier ** p, * q ;
 483   int i, j ;
 484
 485   matrix = (PipMatrix *)malloc(sizeof(PipMatrix)) ;
 486   if (matrix == NULL)
 487   { fprintf(stderr, "Memory Overflow.\n") ;
 488     exit(1) ;
 489   }
 490   matrix->NbRows = NbRows ;
 491   matrix->NbColumns = NbColumns ;
 492   matrix->p_Init_size = NbRows * NbColumns ;
 493   if (NbRows == 0)
 494   { matrix->p = NULL ;
 495     matrix->p_Init = NULL ;
 496   }
 497   else
 498   { if (NbColumns == 0)
 499     { matrix->p = NULL ;
 500       matrix->p_Init = NULL ;
 501     }
 502     else
 503     { p = (Entier **)malloc(NbRows*sizeof(Entier *)) ;
 504       if (p == NULL)
 505       { fprintf(stderr, "Memory Overflow.\n") ;
 506         exit(1) ;
 507       }
 508       q = (Entier *)malloc(NbRows * NbColumns * sizeof(Entier)) ;
 509       if (q == NULL)
 510       { fprintf(stderr, "Memory Overflow.\n") ;
 511         exit(1) ;
 512       }
 513       matrix->p = p ;
 514       matrix->p_Init = q ;
 515       for (i=0;i<NbRows;i++)
 516       { *p++ = q ;
 517         for (j=0;j<NbColumns;j++)
 518         #if defined(LINEAR_VALUE_IS_MP)
 519         mpz_init_set_si(*(q+j),0) ;
 520         #else
 521         *(q+j) = 0 ;
 522         #endif
 523         q += NbColumns ;
 524       }
 525     }
 526   }
 527   return matrix ;
 528 }
 529
 530
 531 /* Fonction pip_matrix_read :
 532  * Adaptation de Matrix_Read de la polylib. Cette fonction lit les donnees
 533  * d'une matrice dans un fichier 'foo' et retourne un pointeur vers une
 534  * structure ou elle a copie les informations de cette matrice. Les donnees
 535  * doivent se presenter tq :
 536  * - des lignes de commentaires commencant par # (optionnelles),
 537  * - le nombre de lignes suivit du nombre de colonnes de la matrice, puis
 538  *   eventuellement d'un commentaire sur une meme ligne,
 539  * - des lignes de la matrice, chaque ligne devant etre sur sa propre ligne de
 540  *   texte et eventuellement suivies d'un commentaire.
 541  * Premiere version : Ced. 18 octobre 2001.
 542  * 24 octobre 2002 : premiere version MP, attention, uniquement capable de
 543  *                   lire des long long pour l'instant. On utilise pas
 544  *                   mpz_inp_str car on lit depuis des char * et non des FILE.
 545  */
 546 PipMatrix * pip_matrix_read(FILE * foo)
 547 { unsigned NbRows, NbColumns ;
 548   int i, j, n ;
 549   #if defined(LINEAR_VALUE_IS_MP)
 550   long long val ;
 551   #endif
 552   char *c, s[1024], str[1024] ;
 553   PipMatrix * matrix ;
 554   Entier * p ;
 555
 556   while (fgets(s,1024,foo) == 0) ;
 557   while ((*s=='#' || *s=='\n') || (sscanf(s," %u %u",&NbRows,&NbColumns)<2))
 558   fgets(s, 1024, foo) ;
 559
 560   matrix = pip_matrix_alloc(NbRows,NbColumns) ;
 561
 562   p = matrix->p_Init ;
 563   for (i=0;i<matrix->NbRows;i++)
 564   { do
 565     { c = fgets(s,1024,foo) ;
 566       while (isspace(*c) && (*c != '\n'))
 567       c++ ;
 568     }
 569     while (c != NULL && (*c == '#' || *c == '\n'));
 570
 571     if (c == NULL)
 572     { fprintf(stderr, "Not enough rows.\n") ;
 573       exit(1) ;
 574     }
 575     for (j=0;j<matrix->NbColumns;j++)
 576     { if (c == NULL || *c == '#' || *c == '\n')
 577       { fprintf(stderr, "Not enough columns.\n") ;
 578         exit(1) ;
 579       }
 580       /* NdCed : Dans le n ca met strlen(str). */
 581       if (sscanf(c,"%s%n",str,&n) == 0)
 582       { fprintf(stderr, "Not enough rows.\n") ;
 583         exit(1) ;
 584       }
 585       #if defined(LINEAR_VALUE_IS_MP)
 586       sscanf(str,"%lld",&val) ;
 587       mpz_set_si(*p++,val) ;
 588       #else
 589       sscanf(str,FORMAT,p++) ;
 590       #endif
 591       c += n ;
 592     }
 593   }
 594   return matrix ;
 595 }
 596
 597
 598 /* initialization of pip */
 599 static int pip_initialized = 0;
 600
 601 void pip_init() {
 602   /* Avoid initializing (and leaking) several times */
 603   if (!pip_initialized) {
 604     #if defined(LINEAR_VALUE_IS_MP)
 605     mpz_init_set_si(UN, 1);
 606     mpz_init_set_si(ZERO, 0);
 607     #else
 608     UN   = VAL_UN ;
 609     ZERO = VAL_ZERO ;
 610     #endif
 611     sol_init() ;
 612     tab_init() ;
 613     pip_initialized = 1;
 614   }
 615 }
 616
 617 void pip_close() {
 618   tab_close();
 619   sol_close();
 620 # if defined(LINEAR_VALUE_IS_MP)
 621   mpz_clear(UN);
 622   mpz_clear(ZERO);
 623 # endif
 624   pip_initialized = 0;
 625 }
 626
 627 /*
 628  * Compute dual variables of equalities from dual variables of
 629  * the corresponding pair of inequalities.
 630  *
 631  * In practice, this means we need to remove one of the two
 632  * dual variables that corrsponding to the inequalities
 633  * and negate the value if it corresponds to the negative
 634  * inequality.
 635  */
 636 static void pip_quast_equalities_dual(PipQuast *solution, PipMatrix *inequnk)
 637 {
 638     PipList **list_p, *list;
 639     int i;
 640
 641     if (!solution)
 642         return;
 643     if (solution->condition) {
 644         pip_quast_equalities_dual(solution->next_then, inequnk);
 645         pip_quast_equalities_dual(solution->next_else, inequnk);
 646     }
 647     if (!solution->list)
 648         return;
 649     if (!solution->next_then)
 650         return;
 651     if (!solution->next_then->list)
 652         return;
 653     list_p = &solution->next_then->list;
 654     for (i = 0; i < inequnk->NbRows; ++i) {
 655         if (entier_zero_p(inequnk->p[i][0])) {
 656             if (entier_notzero_p((*list_p)->vector->the_vector[0])) {
 657                 list_p = &(*list_p)->next;
 658                 list = *list_p;
 659                 *list_p = list->next;
 660                 list->next = NULL;
 661                 pip_list_free(list);
 662             } else {
 663                 list = *list_p;
 664                 *list_p = list->next;
 665                 list->next = NULL;
 666                 pip_list_free(list);
 667                 entier_oppose((*list_p)->vector->the_vector[0],
 668                               (*list_p)->vector->the_vector[0]);
 669                 list_p = &(*list_p)->next;
 670             }
 671         } else
 672             list_p = &(*list_p)->next;
 673     }
 674 }
 675
 676
 677 /******************************************************************************
 678  *                           Fonction de resolution                           *
 679  ******************************************************************************/
 680
 681
 682 /* Fonction pip_solve :
 683  * Cette fonction fait un appel a Pip pour la resolution d'un probleme. Le
 684  * probleme est fourni dans les arguments. Deux matrices de la forme de celles
 685  * utilisees dans la Polylib forment les systemes d'equations/inequations :
 686  * un pour les inconnues, l'autre pour les parametres. Bg est le 'bignum'.
 687  * Le dernier argument contient les options guidant le comportement de PIP.
 688  * Cette fonction retourne la solution sous la forme d'un arbre de structures
 689  * PipQuast.
 690  * 30 juillet 2001 : Premiere version, Ced.
 691  * 18 octobre 2001 : suppression de l'argument Np, le nombre de parametres. Il
 692  *                   est a present deduit de ineqpar. Si ineqpar vaut NULL,
 693  *                   c'est que Np vaut 0. S'il y a des parametres mais pas de
 694  *                   contraintes dessus, ineqpar sera une matrice de 0 lignes
 695  *                   mais du bon nombre de colonnes (Np + 2).
 696  * 27 février 2003 : Verbose est maintenant gradué.
 697  *                  -1 -> silence absolu
 698  *                   0 -> silence relatif
 699  *                   1 -> information sur les coupures dans le cas ou on
 700  *                        cherche une solution entière.
 701  *                   2 -> information sur les pivots et les déterminants
 702  *                   3 -> information sur les tableaux.
 703  *                         Chaque option inclut les précédentes. [paf]
 704  * 15 mars 2003    : passage a une structure d'options.
 705  */
 706 PipQuast * pip_solve(inequnk, ineqpar, Bg, options)
 707 PipMatrix * inequnk, * ineqpar ;
 708 int Bg ;
 709 PipOptions * options ;
 710 { Tableau * ineq, * context, * ctxt ;
 711   int i, Np, Nn, Nl, Nm, p, q, xq, non_vide, Shift = 0, Urs_parms = 0;
 712   char * g ;
 713   struct high_water_mark hq ;
 714   Entier D ;
 715   PipQuast * solution ;
 716   int   sol_flags = 0;
 717
 718   pip_init() ;
 719
 720   /* initialisations diverses :
 721    * - la valeur de Verbose et Deepest_cut sont placees dans leurs variables
 722    *   globales. Dans le cas ou on doit etre en mode verbose, on ouvre le
 723    *   fichier dans lequel ecrire les tracages. Si la variable d'environnement
 724    *   DEBUG est placee, on ecrira dans le nom de fichier correspondant, sinon,
 725    *   dans un nouveau fichier de nom genere par mkstemp,
 726    * - limit est mis au 0 long int (sa valeur par defaut dans Pip original),
 727    * - on lance les initialisations pour tab et sol (autres mises en place
 728    *   de variables globales).
 729    */
 730   verbose = options->Verbose ;
 731   deepest_cut = options->Deepest_cut ;
 732   if (verbose > 0)
 733   { g = getenv("DEBUG") ;
 734     if(g && *g)
 735     { dump = fopen(g, "w") ;
 736       if(dump == NULL)
 737       { fprintf(stderr,"%s unaccessible\n",g) ;
 738         verbose = 0 ;
 739       }
 740     }
 741     else
 742     { mkstemp(dump_name) ;
 743       dump = fopen(dump_name, "w") ;
 744     }
 745   }
 746   #if defined(LINEAR_VALUE_IS_MP)
 747   limit = 0LL ;
 748   #else
 749   limit = ZERO ;
 750   #endif
 751
 752   /* Si inequnk est NULL, la solution est automatiquement void (NULL). */
 753   if (inequnk != NULL)
 754   { /* Np vaut 0 si ineqpar vaut NULL, ineqpar->NbColumns - 2 sinon (on a -1
 755      * pour la constante et -1 pour le marqueur d'egalite/inegalite = -2).
 756      */
 757     Np = (ineqpar == NULL) ? 0 : ineqpar->NbColumns - 2 ;
 758     /* Calcul du nombre d'inconnues du probleme. Comme les matrices d'entree
 759      * sont sous la forme de la polylib.
 760      */
 761     Nn = inequnk->NbColumns - Np - 2 ;
 762     /* Calcul du nombre d'inequations du probleme. Le format de matrice de la
 763      * polylib permet les egalites, on doit donc les compter double quand il
 764      * y en a.
 765      */
 766     Nl = inequnk->NbRows ;
 767     for (i=0;i<inequnk->NbRows;i++)
 768     #if defined(LINEAR_VALUE_IS_MP)
 769     if (mpz_sgn(**(inequnk->p + i)) == 0)
 770     #else
 771     if (**(inequnk->p + i) == 0)
 772     #endif
 773     Nl ++ ;
 774
 775     /* On prend les 'marques' de debut de traitement. */
 776     cross_product = 0 ;
 777     hq = tab_hwm() ;
 778     xq = p = sol_hwm();
 779
 780     if (options->Maximize) {
 781       sol_flags |= SOL_MAX;
 782       Shift = 1;
 783     } else if (options->Urs_unknowns) {
 784       sol_flags |= SOL_SHIFT;
 785       Shift = -1;
 786     }
 787
 788     if (options->Urs_parms) {
 789       Urs_parms = Np - (Bg >= 0);
 790       Np += Urs_parms;
 791     }
 792
 793     /* Si un maximum est demande, mais sans bignum, on crée le bignum. */
 794     if (options->Maximize || options->Urs_unknowns) {
 795       if (Bg < 0) {
 796         Bg = inequnk->NbColumns - 1 ; /* On choisit sa place. */
 797         Np ++ ;                       /* On le compte comme parametre. */
 798         sol_flags |= SOL_REMOVE;      /* On le supprime apres. */
 799       }
 800     }
 801
 802     /* On s'assure d'abord que le systeme pour le contexte n'est pas vide
 803      * avant de commencer le traitement. Si c'est le cas, la solution est
 804      * void (NULL).
 805      */
 806     if (ineqpar != NULL)
 807     { /* Calcul du nombre d'inequations sur les parametres. Le format de
 808        * matrice de la polylib permet les egalites, on doit donc les compter
 809        * double quand il y en a.
 810        */
 811       Nm = ineqpar->NbRows ;
 812       for (i=0;i<ineqpar->NbRows;i++)
 813       #if defined(LINEAR_VALUE_IS_MP)
 814       if (mpz_sgn(**(ineqpar->p + i)) == 0)
 815       #else
 816       if (**(ineqpar->p + i) == 0)
 817       #endif
 818       Nm ++ ;
 819
 820       context = tab_Matrix2Tableau(ineqpar,Nm,Np,0, Shift,Bg-Nn, Urs_parms);
 821       if (options->Nq)
 822         tab_simplify(context, Np);
 823
 824       if (Nm)
 825       { /* Traduction du format de matrice de la polylib vers celui de
 826          * traitement de Pip. Puis traitement proprement dit.
 827          */
 828         ctxt = expanser(context, Np, Nm, Np+1, Np, 0, 0) ;
 829         traiter(ctxt, NULL, Np, 0, Nm, 0, -1, TRAITER_INT);
 830         non_vide = is_not_Nil(p) ;
 831         sol_reset(p) ;
 832       }
 833       else
 834       non_vide = Pip_True ;
 835     }
 836     else
 837     { Nm = 0 ;
 838       ineqpar = pip_matrix_alloc(0, 2);
 839       context = tab_Matrix2Tableau(ineqpar, Nm, Np, 0, Shift, Bg-Nn, Urs_parms);
 840       pip_matrix_free(ineqpar);
 841       non_vide = Pip_True ;
 842     }
 843
 844     if (verbose > 0)
 845     fprintf(dump, "%d %d %d %d %d %d\n",Nn,Np,Nl,Nm,Bg,options->Nq) ;
 846
 847     /* S'il est possible de trouver une solution, on passe au traitement. */
 848     if (non_vide) {
 849       int flags = 0;
 850       ineq = tab_Matrix2Tableau(inequnk,Nl,Nn,Nn, Shift,Bg, Urs_parms);
 851       if (options->Nq)
 852         tab_simplify(ineq, Nn);
 853
 854       compa_count = 0 ;
 855       if (options->Nq)
 856         flags |= TRAITER_INT;
 857       else if (options->Compute_dual) {
 858         flags |= TRAITER_DUAL;
 859         sol_flags |= SOL_DUAL;
 860       }
 861       traiter(ineq, context, Nn, Np, Nl, Nm, Bg, flags);
 862
 863       if (options->Simplify)
 864       sol_simplify(xq) ;
 865       q = sol_hwm() ;
 866       /* On traduit la solution du format de solution de Pip vers un arbre
 867        * de structures de type PipQuast.
 868        */
 869       solution = sol_quast_edit(&xq, NULL, Bg-Nn-1, Urs_parms, sol_flags);
 870       if ((sol_flags & SOL_DUAL) && Nl > inequnk->NbRows)
 871           pip_quast_equalities_dual(solution, inequnk);
 872       sol_reset(p) ;
 873     }
 874     else
 875     return NULL ;
 876     tab_reset(hq) ;
 877   }
 878   else
 879   return NULL ;
 880
 881   return(solution) ;
 882 }