source/piplib.c

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