source/solMP.c

   1 /********************************************************/
   2 /* Part of MultiPrecision PIP port (C) Zbigniew Chamski */
   3 /* and Paul Feautrier.                                  */
   4 /* Based on straight PIP version E.1 by Paul Feautrier  */
   5 /* (<Paul.Feautrier@inria.fr>)                          */
   6 /*                                                      */
   7 /* and a previous port (C) Zbigniew CHAMSKI, 1993.      */
   8 /* (Zbigniew.Chamski@philips.com)                       */
   9 /*                                                      */
  10 /* Copying subject to the terms and conditions of the   */
  11 /* GNU General Public License.                          */
  12 /*                                                      */
  13 /* Send questions, bug reports and fixes to:            */
  14 /*   <Paul.Feautrier@inria.fr>                          */
  15 /********************************************************/
  16
  17 #include <stdio.h>
  18 #include <stdlib.h>
  19 #include <string.h>
  20
  21 #include <piplib/piplib.h>
  22
  23 extern long int cross_product, limit;
  24 extern int verbose;
  25 extern FILE *dump;
  26
  27 struct S
  28     {int flags;
  29      Entier param1, param2;
  30     };
  31
  32 /* In some cases, param1 and param2 are used to store small integers.
  33    These integers are nevertheless converted to bignums in the
  34    interest of simplicity. Since the handling of the solution is
  35    a negligible part of the algorithm, the corresponding overhead is
  36    deemed negligible.
  37 */
  38
  39
  40 #define Free 0
  41 #define Nil  1
  42 #define If   2
  43 #define List 3
  44 #define Form 4
  45 #define New  5
  46 #define Div  6
  47 #define Val  7
  48 #define Error 8
  49
  50 struct S sol_space[SOL_SIZE];
  51 static int sol_free;
  52
  53 /* pgcd has been removed and replaced by mpz_gcd
  54  */
  55
  56 void sol_init(void)
  57 {
  58  sol_free = 0;
  59 }
  60
  61 int sol_hwm()
  62 {
  63  return(sol_free);
  64 }
  65
  66 void sol_reset(p)
  67 int p;
  68 {int i;
  69  if(p<0 || p>=SOL_SIZE)
  70      {fprintf(stdout, "Syserr : sol_reset : Memory allocation error\n");
  71       exit(40);
  72      }
  73  for(i=p; i<sol_free; i++){
  74    mpz_clear(sol_space[i].param1);
  75    mpz_clear(sol_space[i].param2);
  76  }
  77  sol_free = p;
  78 }
  79
  80 struct S *sol_alloc(void)
  81 {struct S *r;
  82  r = sol_space + sol_free;
  83  r->flags = Free;
  84  mpz_init(r->param1);
  85  mpz_init(r->param2);
  86  sol_free++;
  87  if(sol_free >= SOL_SIZE)
  88      {fprintf(stdout, "The solution is too complex! : sol\n");
  89       exit(26);
  90      }
  91      return(r);
  92 }
  93
  94 void sol_nil(void)
  95 {
  96  struct S * r;
  97  r = sol_alloc();
  98  r -> flags = Nil;
  99  if(verbose > 0)
 100    {fprintf(dump, "\nNil");
 101     fflush(dump);
 102   }
 103 }
 104
 105 void sol_error(int c)
 106 {
 107  struct S *r;
 108  r = sol_alloc();
 109  r->flags = Nil;
 110  mpz_set_si(r->param1, c);
 111  if(verbose > 0) {
 112      fprintf(dump, "Erreur %d\n", c);
 113      fflush(dump);
 114      }
 115 }
 116
 117 int is_not_Nil(p)
 118 int p;
 119 {
 120  return(sol_space[p].flags != Nil);
 121 }
 122
 123 void sol_if(void)
 124 {
 125  struct S *r;
 126  r = sol_alloc();
 127  r -> flags = If;
 128  if(verbose > 0) {
 129      fprintf(dump, "\nIf ");
 130      fflush(dump);
 131    }
 132 }
 133
 134 void sol_list(n)
 135 int n;
 136 {struct S * r;
 137  r = sol_alloc();
 138  r->flags = List;
 139  mpz_set_si(r->param1, n);
 140  if(verbose > 0) {
 141      fprintf(dump, "\nList %d ", n);
 142      fflush(dump);
 143 }
 144 }
 145
 146 void sol_forme(l)
 147 int l;
 148 {
 149  struct S *r;
 150  r = sol_alloc();
 151  r -> flags = Form;
 152  mpz_set_ui(r -> param1, l);
 153  if(verbose > 0) {
 154      fprintf(dump, "\nForme %d ", l);
 155      fflush(dump);
 156    }
 157 }
 158
 159 void sol_new(k)
 160 int k;
 161 {
 162  struct S *r;
 163  if(verbose > 0) {
 164      fprintf(dump, "New %d ", k);
 165      fflush(dump);
 166    }
 167  r = sol_alloc();
 168  r -> flags = New;
 169  mpz_set_ui(r -> param1, k);
 170 }
 171
 172 void sol_div()
 173 {
 174  struct S *r;
 175  r = sol_alloc();
 176  r -> flags = Div;
 177  if(verbose > 0) {
 178      fprintf(dump, "Div ");
 179      fflush(dump);
 180    }
 181 }
 182
 183 void sol_val(n, d)
 184 Entier n, d;
 185 {
 186  struct S *r;
 187
 188  r = sol_alloc();
 189  r -> flags = Val;
 190  mpz_set(r -> param1, n);
 191  mpz_set(r -> param2, d);
 192  if(verbose > 0) {
 193    fputs("val(", dump);
 194    mpz_out_str(dump, 10, n);
 195    putc('/', dump);
 196    mpz_out_str(dump, 10, d);
 197    fputs(") ", dump);
 198   }
 199 }
 200
 201 int skip(int);
 202
 203 /* a` partir d'un point de la solution, sauter un objet
 204 bien forme' ainsi qu'un e'ventuel New et pointer sur l'objet
 205 suivant */
 206
 207 int skip_New (int i)
 208 {
 209  if(sol_space[i].flags != New) return i;
 210  i = skip(i+1);      /* sauter le Div */
 211  return i;
 212 }
 213 /* au lancement, i indexe une cellule qui est la te^te d'un objet.
 214    la valeur retourne'e est la te^te de l'objet qui suit. Les
 215    objets de type New sont e'limine's                                */
 216
 217 int skip (int i)
 218 {int n, f;
 219  while((f = sol_space[i].flags) == Free || f == Error) i++;
 220  switch (sol_space[i].flags) {
 221  case Nil : case Val : i++; break;
 222  case New : i = skip_New(i); break;
 223  case If : i = skip(i+1);        /* sauter le pre'dicat */
 224            i = skip(i);          /* sauter le vrai */
 225            i = skip(i); break;   /* sauter le faux */
 226  case List : case Form : n = mpz_get_si(sol_space[i].param1);
 227            i++;
 228            while(n--) i = skip(i);
 229            break;
 230  case Div : i = skip(i+1);       /* sauter la forme */
 231             i = skip(i);         /* sauter le diviseur */
 232             break;
 233  default : fprintf(stdout,
 234               "Syserr : skip : unknown %d\n", sol_space[i].flags);
 235  }
 236  return skip_New(i);
 237 }
 238 /* simplification de la solution : e'limination des constructions
 239    (if p () ()). N'est en service qu'en pre'sence de l'option -z */
 240
 241 void sol_simplify(int i)
 242 {int j, k, l;
 243  if(sol_space[i].flags == If) {
 244      j = skip(i+1);        /* j : de'but de la partie vraie */
 245      k = skip(j);          /* k : de'but de la partie fausse */
 246      sol_simplify(k);
 247      sol_simplify(j);
 248      if(sol_space[j].flags == Nil && sol_space[k].flags == Nil) {
 249          sol_space[i].flags = Nil;
 250          if(k >= sol_free - 1) sol_free = i+1;
 251          else for(l = i+1; l<=k; l++) sol_space[l].flags = Free;
 252        }
 253    }
 254
 255 }
 256 /* e'dition de la solution */
 257
 258 int sol_edit(FILE *foo, int i)
 259 {int j, n;
 260  struct S *p;
 261  Entier N, D, d;
 262  mpz_init(N);
 263  mpz_init(D);
 264  mpz_init(d);
 265  p = sol_space + i;
 266  for(;;) {
 267    if(p->flags == Free) {
 268      p++;
 269      i++;
 270      continue;
 271    }
 272    if(p->flags == New) {
 273      n = mpz_get_si(p->param1);
 274      fprintf(foo, "(newparm %d ", n);
 275      if(verbose>0)fprintf(dump, "(newparm %d ", n);
 276      i = sol_edit(foo, ++i);
 277      p = sol_space +i;
 278      fprintf(foo, ")\n");
 279      if(verbose>0)fprintf(dump, ")\n");
 280      continue;
 281    }
 282    break;
 283  }
 284  switch(p->flags)
 285    {
 286    case Nil : fprintf(foo, "()\n");
 287      if(verbose>0)fprintf(dump, "()\n");
 288      i++; break;
 289    case Error : n = mpz_get_si(p->param1);
 290      fprintf(foo, "Error %d\n", n);
 291      if(verbose>0)fprintf(dump, "Error %d\n", n);
 292      i++; break;
 293    case If  : fprintf(foo, "(if ");
 294      if(verbose>0)fprintf(dump, "(if ");
 295      i = sol_edit(foo, ++i);
 296      i = sol_edit(foo, i);
 297      i = sol_edit(foo, i);
 298      fprintf(foo, ")\n");
 299      if(verbose>0)fprintf(dump, ")\n");
 300      break;
 301    case List: fprintf(foo, "(list ");
 302      if(verbose>0)fprintf(dump, "(list ");
 303      n = mpz_get_si(p->param1);
 304      i++;
 305      while(n--) i = sol_edit(foo, i);
 306      fprintf(foo, ")\n");
 307      if(verbose>0)fprintf(dump, ")\n");
 308      break;
 309    case Form: fprintf(foo, "#[");
 310      if(verbose>0)fprintf(dump, "#[");
 311      n = mpz_get_si(p->param1);
 312      for(j = 0; j<n; j++){
 313        i++; p++;
 314        mpz_set(N, p->param1); mpz_set(D, p->param2);
 315        mpz_gcd(d, N, D);
 316        if(mpz_cmp(d, D) == 0){
 317          putc(' ', foo);
 318          mpz_divexact(N, N, d);
 319          mpz_out_str(foo, 10, N);
 320          if(verbose>0){
 321            putc(' ', dump);
 322            mpz_out_str(dump, 10, N);
 323          }
 324        }
 325        else{
 326          mpz_divexact(N, N, d);
 327          mpz_divexact(D, D, d);
 328          putc(' ', foo);
 329          mpz_out_str(foo, 10, N);
 330          putc('/', foo);
 331          mpz_out_str(foo, 10, D);
 332          if(verbose>0){
 333            putc(' ', dump);
 334            mpz_out_str(dump, 10, N);
 335            putc('/', dump);
 336            mpz_out_str(dump, 10, D);
 337          }
 338        }
 339      }
 340      fprintf(foo, "]\n");
 341      if(verbose>0)
 342        fputs("]\n", dump);
 343      i++;
 344      break;
 345    case Div : fprintf(foo, "(div ");
 346      if(verbose>0)fprintf(dump, "(div ");
 347      i = sol_edit(foo, ++i);
 348      i = sol_edit(foo, i);
 349      fprintf(foo, ")\n");
 350      if(verbose>0)fprintf(dump, ")\n");
 351      break;
 352    case Val : mpz_set(N, p->param1); mpz_set(D, p->param2);
 353      mpz_gcd(d, N, D);
 354      if(mpz_cmp(d, D) == 0){
 355        mpz_divexact(N, N, d);
 356        putc(' ', foo);
 357        mpz_out_str(foo, 10, N);
 358        if(verbose>0){
 359          putc(' ', dump);
 360          mpz_out_str(dump, 10, N);
 361        }
 362      }
 363      else{
 364        mpz_divexact(N, N, d);
 365        mpz_divexact(D, D, d);
 366        putc(' ', foo);
 367        mpz_out_str(foo, 10, N);
 368        fprintf(foo, "/");
 369        mpz_out_str(foo, 10, D);
 370        if(verbose>0){
 371          putc(' ', dump);
 372          mpz_out_str(dump, 10, N);
 373          fprintf(dump, "/");
 374          mpz_out_str(dump, 10, D);
 375        }
 376      }
 377      i++;
 378      break;
 379    default  : fprintf(foo, "Inconnu : sol\n");
 380      if(verbose>0)fprintf(dump, "Inconnu : sol\n");
 381    }
 382  mpz_clear(d);
 383  mpz_clear(D);
 384  mpz_clear(N);
 385  return(i);
 386 }
 387
 388
 389
 390
 391
 392 /* Fonction sol_vector_edit :
 393  * Cette fonction a pour but de placer les informations correspondant
 394  * a un Vector dans la grammaire dans une structure de type PipVector. Elle
 395  * prend en parametre un pointeur vers une case memoire contenant le
 396  * numero de cellule du tableau sol_space a partir de laquelle on doit
 397  * commencer la lecture des informations. Elle retourne un pointeur vers
 398  * une structure de type PipVector contenant les informations de ce Vector.
 399  * Premiere version : Ced. 20 juillet 2001.
 400  * 24 octobre 2002 : Premiere version MP.
 401  */
 402 PipVector * sol_vector_edit(int * i)
 403 { int j, n ;
 404   struct S *p ;
 405   Entier N, D, d ;
 406   PipVector * vector ;
 407
 408   mpz_init(N) ;
 409   mpz_init(D) ;
 410   mpz_init(d) ;
 411
 412   vector = (PipVector *)malloc(sizeof(PipVector)) ;
 413   if (vector == NULL)
 414   { fprintf(stderr, "Memory Overflow.\n") ;
 415     exit(1) ;
 416   }
 417   p = sol_space + (*i) ;
 418   n = mpz_get_si(p->param1) ;
 419   vector->nb_elements = n ;
 420   vector->the_vector = (Entier *)malloc(sizeof(Entier)*n) ;
 421   if (vector->the_vector == NULL)
 422   { fprintf(stderr, "Memory Overflow.\n") ;
 423     exit(1) ;
 424   }
 425   vector->the_deno = (Entier *)malloc(sizeof(Entier)*n) ;
 426   if (vector->the_deno == NULL)
 427   { fprintf(stderr, "Memory Overflow.\n") ;
 428     exit(1) ;
 429   }
 430
 431   for (j=0;j<n;j++)
 432   { (*i)++ ;
 433     p++ ;
 434     mpz_set(N,p->param1) ;
 435     mpz_set(D,p->param2) ;
 436     mpz_gcd(d, N, D);
 437
 438     mpz_init(vector->the_vector[j]) ;
 439     mpz_divexact(vector->the_vector[j],N,d) ;
 440
 441     mpz_init(vector->the_deno[j]) ;
 442     if (mpz_cmp(d, D) == 0)
 443     mpz_set(vector->the_deno[j],UN) ;
 444     else
 445     mpz_divexact(vector->the_deno[j],D,d) ;
 446   }
 447   (*i)++ ;
 448
 449   mpz_clear(d);
 450   mpz_clear(D);
 451   mpz_clear(N);
 452   return(vector) ;
 453 }
 454
 455
 456 /* Fonction sol_newparm_edit :
 457  * Cette fonction a pour but de placer les informations correspondant
 458  * a un Newparm dans la grammaire dans une structure de type PipNewparm. Elle
 459  * prend en parametre un pointeur vers une case memoire contenant le
 460  * numero de cellule du tableau sol_space a partir de laquelle on doit
 461  * commencer la lecture des informations. Elle retourne un pointeur vers
 462  * une structure de type PipNewparm contenant les informations de ce Newparm.
 463  * Premiere version : Ced. 18 octobre 2001.
 464  * 24 octobre 2002 : Premiere version MP.
 465  */
 466 PipNewparm * sol_newparm_edit(int * i)
 467 { struct S * p ;
 468   PipNewparm * newparm, * newparm_new, * newparm_now ;
 469
 470   /* On place p au lieu de lecture. */
 471   p = sol_space + (*i) ;
 472   /* On passe le New et le Div pour aller a Form et lire le VECTOR. */
 473   (*i) += 2 ;
 474
 475   newparm = (PipNewparm *)malloc(sizeof(PipNewparm)) ;
 476   if (newparm == NULL)
 477   { fprintf(stderr, "Memory Overflow.\n") ;
 478     exit(1) ;
 479   }
 480   newparm->vector = sol_vector_edit(i) ;
 481   newparm->rank = mpz_get_si(p->param1) ;
 482   /* On met p a jour pour lire le denominateur (un Val de param2 UN). */
 483   p = sol_space + (*i) ;
 484   mpz_init_set(newparm->deno,p->param1) ;
 485   newparm->next = NULL ;
 486
 487   newparm_now = newparm ;
 488   if (verbose)
 489   { fprintf(dump,"\n(newparm ") ;
 490     fprintf(dump,"%d",newparm->rank) ;
 491     fprintf(dump," (div ") ;
 492     pip_vector_print(dump,newparm->vector) ;
 493     fprintf(dump," ") ;
 494     mpz_out_str(dump,10,newparm->deno) ;
 495     fprintf(dump,"))") ;
 496   }
 497
 498   /* On passe aux elements suivants. */
 499   (*i) ++ ;
 500   p = sol_space + (*i) ;
 501   while (p->flags == New)
 502   { (*i) += 2 ;
 503     newparm_new = (PipNewparm *)malloc(sizeof(PipNewparm)) ;
 504     if (newparm_new == NULL)
 505     { fprintf(stderr, "Memory Overflow.\n") ;
 506       exit(1) ;
 507     }
 508     newparm_new->vector = sol_vector_edit(i) ;
 509     newparm_new->rank = mpz_get_si(p->param1) ;
 510     p = sol_space + (*i) ;
 511     mpz_init_set(newparm_new->deno,p->param1) ;
 512     newparm_new->next = NULL ;
 513
 514     newparm_now->next = newparm_new ;
 515     newparm_now = newparm_now->next ;
 516     if (verbose)
 517     { fprintf(dump,"\n(newparm ") ;
 518       fprintf(dump,"%d",newparm_new->rank) ;
 519       fprintf(dump," (div ") ;
 520       pip_vector_print(dump,newparm_new->vector) ;
 521       fprintf(dump," ") ;
 522       mpz_out_str(dump,10,newparm_new->deno) ;
 523       fprintf(dump,"))") ;
 524     }
 525     (*i) ++ ;
 526     p = sol_space + (*i) ;
 527   }
 528   return(newparm) ;
 529 }
 530
 531
 532 /* Fonction sol_list_edit :
 533  * Cette fonction a pour but de placer les informations correspondant
 534  * a une List dans la grammaire dans une structure de type PipList. Elle
 535  * prend en parametre un pointeur vers une case memoire contenant le
 536  * numero de cellule du tableau sol_space a partir de laquelle on doit
 537  * commencer la lecture des informations. Elle retourne un pointeur vers
 538  * une structure de type PipList contenant les informations de cette List.
 539  * Premiere version : Ced. 18 octobre 2001.
 540  */
 541 PipList * sol_list_edit(int * i, int nb_elements)
 542 { PipList * list, * list_new, * list_now ;
 543
 544   /* Pour le premier element. */
 545   list = (PipList *)malloc(sizeof(PipList)) ;
 546   if (list == NULL)
 547   { fprintf(stderr, "Memory Overflow.\n") ;
 548     exit(1) ;
 549   }
 550   list->vector = sol_vector_edit(i) ;
 551   list->next = NULL ;
 552
 553   list_now = list ;
 554   if (verbose)
 555   { fprintf(dump,"\n(list ") ;
 556     pip_vector_print(dump,list->vector) ;
 557   }
 558   nb_elements-- ;
 559
 560   /* Pour les elements suivants. */
 561   while (nb_elements--)
 562   { list_new = (PipList *)malloc(sizeof(PipList)) ;
 563     if (list_new == NULL)
 564     { fprintf(stderr, "Memory Overflow.\n") ;
 565       exit(1) ;
 566     }
 567     list_new->vector = sol_vector_edit(i) ;
 568     list_new->next = NULL ;
 569
 570     if (verbose)
 571     { fprintf(dump,"\n") ;
 572       pip_vector_print(dump,list_new->vector) ;
 573     }
 574     list_now->next = list_new ;
 575     list_now = list_now->next ;
 576   }
 577   if (verbose)
 578   fprintf(dump,"\n)") ;
 579
 580   return(list) ;
 581 }
 582
 583
 584 /* Fonction sol_quast_edit :
 585  * Cette fonction a pour but de placer les informations de la solution
 586  * (qui sont contenues dans le tableau sol_space) dans une structure de
 587  * type PipQuast en vue d'une utilisation directe de la solution par une
 588  * application exterieure. Elle prend en parametre un pointeur vers une
 589  * case memoire contenant le numero de cellule du tableau sol_space
 590  * a partir de laquelle on doit commencer la lecture des informations. Elle
 591  * recoit aussi l'adresse du PipQuast qui l'a appelle (pour le champ father).
 592  * Elle retourne un pointeur vers une structure de type PipQuast qui
 593  * contient toutes les informations sur la solution (sous forme d'arbre).
 594  * Remarques : cette fonction lit les informations comme elles doivent
 595  * se presenter a la fin du traitement. Elle respecte scrupuleusement
 596  * la grammaire attendue et n'accepte de passer des cellules a Free
 597  * qu'entre une des trois grandes formes (if, list ou suite de newparm).
 598  * 20 juillet 2001 : Premiere version, Ced.
 599  * 31 juillet 2001 : Ajout du traitement de l'option verbose = code*2 :0(
 600  * 18 octobre 2001 : Grands changements dus a l'eclatement de la structure
 601  *                   PipVector en PipVector, PipNewparm et PipList, et
 602  *                   eclatement de la fonction avec sol_newparm_edit et
 603  *                   sol_list_edit.
 604  * 24 octobre 2002 : Premiere version MP.
 605  */
 606 PipQuast * sol_quast_edit(int * i, PipQuast * father)
 607 { int nb_elements ;
 608   struct S * p ;
 609   PipQuast * solution ;
 610   PipList * list_new, * list_now ;
 611   PipNewparm * newparm_new, * newparm_now ;
 612
 613   /* On place p au lieu de lecture. */
 614   p = sol_space + (*i) ;
 615   /* En cas d'utilisation de l'option de simplification, une plage de
 616    * structures S peut avoir les flags a Free. On doit alors les passer.
 617    */
 618   while (p->flags == Free)
 619   { p ++ ;
 620     (*i) ++ ;
 621   }
 622
 623   solution = (PipQuast *)malloc(sizeof(PipQuast)) ;
 624   if (solution == NULL)
 625   { fprintf(stderr, "Memory Overflow.\n") ;
 626     exit(1) ;
 627   }
 628   solution->newparm = NULL ;
 629   solution->list = NULL ;
 630   solution->condition = NULL ;
 631   solution->next_then = NULL ;
 632   solution->next_else = NULL ;
 633   solution->father = father ;
 634
 635   /* On peut commencer par une chaine de nouveaux parametres... */
 636   if (p->flags == New)
 637   { solution->newparm = sol_newparm_edit(i) ;
 638     p = sol_space + (*i) ;
 639   }
 640
 641   /* ...ensuite soit par une liste (vide ou non) soit par un if. */
 642   (*i)++ ; /* Factorise de List, Nil et If. */
 643   switch (p->flags)
 644   { case List : if ((nb_elements = mpz_get_si(p->param1)) != 0)
 645                 solution->list = sol_list_edit(i,nb_elements) ;
 646                 break ;
 647     case Nil  : if (verbose)
 648                 fprintf(dump,"\n()") ;
 649                 break ;
 650     case If   : solution->condition = sol_vector_edit(i) ;
 651                 if (verbose)
 652                 { fprintf(dump,"\n(if ") ;
 653                   pip_vector_print(dump,solution->condition) ;
 654                 }
 655                 solution->next_then = sol_quast_edit(i,solution) ;
 656                 solution->next_else = sol_quast_edit(i,solution) ;
 657                 if (verbose)
 658                 fprintf(dump,"\n)") ;
 659                 break ;
 660     default   : fprintf(stderr,"\nAie !!! Flag %d inattendu.\n",p->flags) ;
 661                 if (verbose)
 662                 fprintf(dump,"\nAie !!! Flag %d inattendu.\n",p->flags) ;
 663                 exit(1) ;
 664   }
 665
 666   return(solution) ;
 667 }
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692