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