source/traiterMP.c

   1 /********************************************************/
   2 /* Part of MultiPrecision PIP port by Zbigniew Chamski  */
   3 /* and Paul Feautrier.                                  */
   4 /* Based on straight PIP E.1 version 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
  20 #include <piplib/piplib.h>
  21
  22 #define max(x,y) ((x) > (y)? (x) : (y))
  23
  24 extern long int cross_product, limit;
  25 extern int verbose;
  26 extern FILE *dump;
  27 extern int profondeur;
  28 extern float clock;
  29 extern int compa_count;
  30
  31 /* The function llog is replaced by mpz_sizeinbase. */
  32
  33 int chercher(Tableau *p, int masque, int n)
  34 {int i;
  35  for(i = 0; i<n; i++)
  36      if(p->row[i].flags & masque) break;
  37  return(i);
  38 }
  39
  40 /* il est convenu que traiter ne doit modifier ni le tableau, ni le contexte;
  41    le tableau peut grandir en cas de coupure (+1 en hauteur et +1 en largeur
  42    si nparm != 0) et en cas de partage (+1 en hauteur)
  43    (seulement si nparm != 0).
  44    le contexte peut grandir en cas de coupure (+2 en hauteur et +1 en largeur)
  45    (seulement si nparm !=0) et en cas de partage (+1 en hauteur)(nparm !=0).
  46    On estime le nombre de coupures a llog(D) et le nombre de partages a
  47    ni.
  48 */
  49
  50 Tableau *expanser(Tableau *tp, int virt, int reel, int ncol,
  51                                int off, int dh, int dw)
  52 {
  53  int i, j, ff;
  54  char *q; Entier *pq;
  55  Entier *pp, *qq;
  56  Tableau *rp;
  57  if(tp == NULL) return(NULL);
  58  rp = tab_alloc(reel+dh, ncol+dw, virt);
  59
  60  mpz_set(rp->determinant, tp->determinant);
  61  pq = (Entier *) & (rp->row[virt+reel+dh]);
  62  for(i = off; i<virt + reel; i++)
  63      {ff = Flag(rp, i) = Flag(tp, i-off);
  64       mpz_set(Denom(rp, i), Denom(tp, i-off));
  65       if(ff & Unit) rp->row[i].objet.unit = tp->row[i-off].objet.unit;
  66       else {
  67           rp->row[i].objet.val = pq;
  68           pq +=(ncol + dw);
  69           pp = tp->row[i-off].objet.val;
  70           qq = rp->row[i].objet.val;
  71           for(j = 0; j<ncol; j++) mpz_set(*qq++, *pp++);
  72           }
  73       }
  74  return(rp);
  75 }
  76
  77 int exam_coef(Tableau *tp, int nvar, int ncol, int bigparm)
  78 {int i, j;
  79  int ff, fff;
  80  int x;
  81  Entier *p;
  82  for(i = 0; i<tp->height; i++)
  83      {ff = Flag(tp,i);
  84       if(ff == 0) break;
  85       if(ff == Unknown)
  86           {if(bigparm >= 0){
  87             x = mpz_sgn(Index(tp,i, bigparm));
  88             if(x<0){
  89               Flag(tp, i) = Minus;
  90               return(i);
  91             }
  92             else if(x>0)
  93               Flag(tp, i) = Plus;
  94             continue;
  95           }
  96           ff = Zero;
  97           p = &(tp->row[i].objet.val[nvar+1]);
  98           for(j = nvar+1; j<ncol; j++){
  99             x = mpz_sgn(*p); p++;
 100             if(x<0) fff = Minus;
 101             else if (x>0) fff = Plus;
 102             else fff = Zero;
 103             if(fff != Zero && fff != ff)
 104               if(ff == Zero) ff = fff;
 105               else {ff = Unknown;
 106               break;
 107               }
 108           }
 109 /* bug de'tecte' par [paf], 16/2/93 !
 110    Si tous les coefficients des parame`tres sont ne'gatifs
 111    et si le terme constant est nul, le signe est inconnu!!
 112    On traite donc spe'cialement le terme constant. */
 113           x = mpz_sgn(Index(tp, i, nvar));
 114           if(x<0) fff = Minus;
 115           else if(x>0) fff = Plus;
 116           else fff = Zero;
 117 /* ici on a le signe du terme constant */
 118           switch(ff){
 119 /* le signe est inconnu si les coefficients sont positifs et
 120    le terme constant ne'gatif */
 121           case Plus: if(fff == Minus) ff = Unknown; break;
 122 /* si les coefficients sont tous nuls, le signe est celui
 123    du terme constant */
 124           case Zero: ff = fff; break;
 125 /* le signe est inconnu si les coefficients sont ne'gatifs,
 126    sauf si le terme constant est ‚galement n‚gatif. */
 127           case Minus: if(fff != Minus) ff = Unknown; break;
 128 /* enfin, il n'y a rien a` dire si le signe des coefficients est inconnu */
 129           }
 130           Flag(tp, i) = ff;
 131           if(ff == Minus) return(i);
 132           }
 133      }
 134  return(i);
 135 }
 136
 137 void compa_test(Tableau *tp, Tableau *context,
 138                 int ni, int nvar, int nparm, int nc)
 139 {
 140  Entier discr[MAXPARM];
 141  int i, j;
 142  int ff;
 143  int cPlus, cMinus, isCritic;
 144  int verbold;
 145  Tableau *tPlus, *tMinus;
 146  int p;
 147  struct high_water_mark q;
 148  if(nparm == 0) return;
 149  if(nparm >= MAXPARM) {
 150      fprintf(stdout, "Too much parameters : %d\n", nparm);
 151      exit(1);
 152      }
 153  q = tab_hwm();
 154  for(i=0; i<=nparm; i++)
 155    mpz_init(discr[i]);
 156
 157  for(i = 0; i<ni + nvar; i++)
 158      {ff = Flag(tp,i);
 159       if(ff & (Critic | Unknown))
 160           {isCritic = True;
 161            for(j = 0; j<nvar; j++) if(mpz_sgn(Index(tp, i, j)) > 0)
 162                  {isCritic = False;
 163                   break;
 164                  }
 165            compa_count++;
 166            for(j = 0; j<nparm; j++)mpz_set(discr[j], Index(tp, i, j+nvar+1));
 167            mpz_set(discr[nparm], Index(tp, i, nvar));
 168            mpz_sub_ui(discr[nparm], discr[nparm], (isCritic ? 0 : 1));
 169            tPlus = expanser(context, nparm, nc, nparm+1, nparm, 1, 0);
 170            Flag(tPlus, nparm+nc) = Unknown;
 171            for(j = 0;j<=nparm; j++)mpz_set(Index(tPlus, nparm+nc, j),discr[j]);
 172            mpz_set(Denom(tPlus, nparm+nc), UN);
 173            p = sol_hwm();
 174            traiter(tPlus, NULL, True, nparm, 0, nc+1, 0, -1);
 175            cPlus = is_not_Nil(p);
 176            if(verbose>0){
 177              fprintf(dump, "\nThe positive case has been found ");
 178              fprintf(dump, cPlus? "possible\n": "impossible\n");
 179              fflush(dump);
 180            }
 181
 182            sol_reset(p);
 183            for(j = 0; j<nparm+1; j++) mpz_neg(discr[j], discr[j]);
 184            mpz_sub_ui(discr[nparm], discr[nparm], (isCritic ? 1 : 2));
 185            tMinus = expanser(context, nparm, nc, nparm+1, nparm, 1, 0);
 186            Flag(tMinus, nparm+nc) = Unknown;
 187            for(j=0;j<=nparm; j++)mpz_set(Index(tMinus, nparm+nc, j),discr[j]);
 188            mpz_set(Denom(tMinus, nparm+nc), UN);
 189            traiter(tMinus, NULL, True, nparm, 0, nc+1, 0, -1);
 190            cMinus = is_not_Nil(p);
 191            if(verbose>0){
 192              fprintf(dump, "\nThe negative case has been found ");
 193              fprintf(dump, cMinus? "possible\n": "impossible\n");
 194              fflush(dump);
 195            }
 196            sol_reset(p);
 197            if (cPlus && cMinus) {
 198                Flag(tp,i) = isCritic ? Critic : Unknown;
 199              }
 200            else if (cMinus)
 201               {Flag(tp,i) = Minus;
 202                break;
 203               }
 204            else {
 205              Flag(tp,i) = cPlus ? Plus : Zero;
 206            }
 207           }
 208      }
 209  tab_reset(q);
 210
 211  for(i=0; i<=nparm; i++)
 212    mpz_clear(discr[i]);
 213
 214  return;
 215 }
 216
 217 Entier *valeur(Tableau *tp, int i, int j)
 218 {
 219  if(Flag(tp, i) & Unit)
 220      return(tp->row[i].objet.unit == j ? &Denom(tp,i) : &ZERO);
 221  else return(&Index(tp, i, j));
 222 }
 223
 224 void solution(Tableau *tp, int nvar, int nparm)
 225 {int i, j;
 226  int ncol = nvar + nparm + 1;
 227  Entier *denom;
 228
 229  sol_list(nvar);
 230  for(i = 0; i<nvar; i++)
 231      {sol_forme(nparm+1);
 232       for(j = nvar+1; j<ncol; j++)
 233          sol_val(*valeur(tp, i, j), Denom(tp,i));
 234       sol_val(*valeur(tp, i, nvar), Denom(tp,i));
 235      }
 236 }
 237
 238 int choisir_piv(Tableau *tp, int pivi, int nvar, int nligne)
 239 {
 240  int j, k;
 241  Entier pivot, foo, x, y;
 242  int sgn_x;
 243  int pivj = -1;
 244
 245  mpz_init(pivot); mpz_init(foo); mpz_init(x); mpz_init(y);
 246
 247  for(j = 0; j<nvar; j++){
 248    mpz_set(foo, Index(tp, pivi, j));
 249    if(mpz_sgn(foo) <= 0) continue;
 250    if(pivj < 0){
 251      pivj = j;
 252      mpz_set(pivot, foo);
 253      continue;
 254    }
 255    for(k = 0; k<nligne; k++){
 256      mpz_mul(x, pivot, *valeur(tp, k, j));
 257      mpz_mul(y, *valeur(tp, k, pivj), foo);
 258      mpz_sub(x, x, y);
 259      cross_product++;
 260      sgn_x = mpz_sgn(x);
 261      if(sgn_x) break;
 262    }
 263    if(sgn_x < 0){
 264      pivj = j;
 265      mpz_set(pivot, foo);
 266    }
 267  }
 268  mpz_clear(pivot); mpz_clear(foo); mpz_clear(x); mpz_clear(y);
 269  return(pivj);
 270 }
 271
 272 int pivoter(Tableau *tp, int pivi, int nvar,
 273                int nparm, int ni, int iq)
 274 {int pivj;
 275  int ncol = nvar + nparm + 1;
 276  int nligne = nvar + ni;
 277  int i, j, k;
 278  Entier x, y, d, gcd, u, dpiv;
 279  int ff, fff;
 280  Entier pivot, foo, z;
 281  Entier ppivot, dppiv;
 282  Entier new[MAXCOL], *p, *q;
 283  Entier lpiv;
 284  int sgn_x;
 285
 286
 287  if(ncol >= MAXCOL) {
 288    fprintf(stdout, "Too much variables\n");
 289    exit(1);
 290  }
 291  if(0 > pivi || pivi >= nligne || Flag(tp, pivi) == Unit) {
 292    fprintf(stdout, "Syserr : pivoter : wrong pivot row\n");
 293    exit(1);
 294  }
 295
 296  pivj = choisir_piv(tp, pivi, nvar, nligne);
 297  if(pivj < 0) return(-1);
 298  if(pivj >= nvar) {
 299    fprintf(stdout, "Syserr : pivoter : wrong pivot\n");
 300    exit(1);
 301  }
 302
 303  mpz_init(x); mpz_init(y); mpz_init(d);
 304  mpz_init(gcd); mpz_init(u); mpz_init(dpiv);
 305  mpz_init(lpiv); mpz_init(pivot); mpz_init(foo);
 306  mpz_init(z); mpz_init(ppivot); mpz_init(dppiv);
 307
 308  for(i=0; i<ncol; i++)
 309    mpz_init(new[i]);
 310
 311  mpz_set(pivot, Index(tp, pivi, pivj));
 312  mpz_set(dpiv, Denom(tp, pivi));
 313  mpz_gcd(d, pivot, dpiv);
 314  mpz_divexact(ppivot, pivot, d);
 315  mpz_divexact(dppiv, dpiv, d);
 316  if(verbose>0){
 317    fprintf(dump, "Pivot ");
 318    mpz_out_str(dump, 10, ppivot);
 319    putc('/', dump);
 320    mpz_out_str(dump, 10, dppiv);
 321    putc('\n', dump);
 322    fprintf(dump, "%d x %d\n", pivi, pivj);
 323  }
 324  mpz_fdiv_qr(x, y, tp->determinant, dppiv);
 325
 326  if(mpz_sgn(y) != 0){
 327    fprintf(stdout, "Computation error\n");
 328    if(verbose>0) fflush(dump);
 329    exit(1);
 330  }
 331  mpz_mul(tp->determinant, x, ppivot);
 332
 333  if(verbose>0){
 334    fprintf(dump, "determinant ");
 335    mpz_out_str(dump, 10, tp->determinant);
 336    fprintf(dump, "\n");
 337  }
 338
 339
 340  for(j = 0; j<ncol; j++)
 341    if(j==pivj)
 342      mpz_set(new[j], dpiv);
 343    else
 344      mpz_neg(new[j], Index(tp, pivi, j));
 345  for(k = 0; k<nligne; k++){
 346    if(Flag(tp,k) & Unit)continue;
 347    if(k == pivi)continue;
 348    /* foo = Index(tp, k, pivj) */
 349    mpz_set(foo, Index(tp, k, pivj));
 350    /* d = gcd(pivot, foo); */
 351    mpz_gcd(d, pivot, foo);
 352    /* lpiv = pivot/d; */
 353    mpz_divexact(lpiv, pivot, d);
 354    /* foo /= d; */
 355    mpz_divexact(foo, foo, d);
 356    /* d = Denom(tp,k); */
 357    mpz_set(d, Denom(tp,k));
 358    mpz_mul(gcd, lpiv, d);
 359    mpz_set(Denom(tp, k), gcd);
 360    p = tp->row[k].objet.val;
 361    q = tp->row[pivi].objet.val;
 362    for(j = 0; j<ncol; j++){
 363      if(j == pivj)
 364        mpz_mul(z, dpiv, foo);
 365      else {
 366        mpz_mul(z, *p, lpiv);
 367        mpz_mul(y, *q, foo);
 368        mpz_sub(z, z, y);
 369      }
 370      q++;
 371      cross_product++;
 372      mpz_set(*p, z);
 373      p++;
 374      if(mpz_cmp_ui(gcd, 1) != 0)
 375        mpz_gcd(gcd, gcd, z);
 376    }
 377    if(mpz_cmp_ui(gcd, 1) != 0){
 378      p = tp->row[k].objet.val;
 379      for(j = 0; j<ncol; j++){
 380        mpz_divexact(*p, *p, gcd);
 381        p++;
 382      }
 383    }
 384    mpz_divexact(Denom(tp,k), Denom(tp,k), gcd);
 385  }
 386  p = tp->row[pivi].objet.val;
 387  for(k = 0; k<nligne; k++)
 388    if((Flag(tp, k) & Unit) && tp->row[k].objet.unit == pivj) break;
 389  Flag(tp, k) = Plus;
 390  tp->row[k].objet.val = p;
 391  for(j = 0; j<ncol; j++)
 392    /* *p++ = new[j] */
 393    mpz_set(*p++, new[j]);
 394  mpz_set(Denom(tp, k), pivot);
 395  Flag(tp, pivi) = Unit | Zero;
 396  mpz_set(Denom(tp, pivi), UN);
 397  tp->row[pivi].objet.unit = pivj;
 398
 399  for(k = 0; k<nligne; k++){
 400    ff = Flag(tp, k);
 401    if(ff & Unit) continue;
 402    sgn_x = mpz_sgn(Index(tp, k, pivj));
 403    if(sgn_x < 0) fff = Minus;
 404    else if(sgn_x == 0) fff = Zero;
 405    else fff = Plus;
 406    if(fff != Zero && fff != ff)
 407      if(ff == Zero) ff = (fff == Minus ? Unknown : fff);
 408      else ff = Unknown;
 409    Flag(tp, k) = ff;
 410  }
 411
 412  if(verbose>0){
 413    fprintf(dump, "just pivoted\n");
 414    tab_display(tp, dump);
 415  }
 416
 417  mpz_clear(x); mpz_clear(y); mpz_clear(d); mpz_clear(gcd);
 418  mpz_clear(u); mpz_clear(dpiv); mpz_clear(lpiv);
 419  mpz_clear(pivot); mpz_clear(foo); mpz_clear(z);
 420  mpz_clear(ppivot); mpz_clear(dppiv);
 421
 422  for(i=0; i<ncol; i++)
 423    mpz_clear(new[i]);
 424
 425  return(0);
 426 }
 427
 428 /* dans cette version, "traiter" modifie ineq; par contre
 429    le contexte est immediatement recopie' */
 430
 431 void traiter(tp, ctxt, iq, nvar, nparm, ni, nc, bigparm)
 432 Tableau *tp, *ctxt;
 433 int iq, nvar, nparm, ni, nc, bigparm;
 434
 435 {
 436  int j;
 437  int pivi, nligne, ncol;
 438  struct high_water_mark x;
 439  Tableau *context;
 440  int dch, dcw;
 441  double s, t, d, smax;
 442  int i;
 443  int sgn_x;
 444  struct L temp;
 445  Entier discr[MAXPARM];
 446  for(i=0; i<=MAXPARM; i++)
 447    mpz_init(discr[i]);
 448
 449  dcw = mpz_sizeinbase(tp->determinant, 2);
 450  dch = 2 * dcw + 1;
 451  x = tab_hwm();
 452  nligne = nvar+ni;
 453
 454  context = expanser(ctxt, 0, nc, nparm+1, 0, dch, dcw);
 455 /*
 456  sort the rows in increasing order of the largest coefficient
 457 */
 458  smax = 0.;
 459
 460  for(i=nvar; i<nligne; i++){
 461    if(Flag(tp,i) & Unit) continue;
 462    s = 0.;
 463    d = mpz_get_d(Denom(tp,i));
 464    for(j=0; j<nvar; j++){
 465      t = mpz_get_d(Index(tp,i,j))/d;
 466      s = max(s, abs(t));
 467    }
 468    tp->row[i].size = s;
 469    smax = max(s, smax);
 470  }
 471
 472  for(i=nvar; i<nligne; i++){
 473    if(Flag(tp,i) & Unit) continue;
 474    s = smax;
 475    pivi = i;
 476    for(j=i; j<nligne; j++){
 477      if(Flag(tp,j) & Unit) continue;
 478      if(tp->row[j].size < s){
 479        s = tp->row[i].size;
 480        pivi = j;
 481      }
 482    }
 483    if(pivi != i) {
 484      temp = tp->row[pivi];
 485      tp->row[pivi] = tp->row[i];
 486      tp->row[i]=temp;
 487    }
 488  }
 489
 490
 491  if(verbose>0){
 492    fprintf(dump, "just sorted\n");
 493    tab_display(tp, dump);
 494    fflush(dump);
 495  }
 496
 497  for(;;) {
 498    nligne = nvar+ni; ncol = nvar+nparm+1;
 499    if(nligne > tp->height || ncol > tp->width) {
 500      fprintf(stdout, "Syserr : traiter : tableau too small\n");
 501      exit(1);
 502    }
 503    pivi = chercher(tp, Minus, nligne);
 504    if(pivi < nligne) goto pirouette;           /* There is a negative row   */
 505
 506    pivi = exam_coef(tp, nvar, ncol, bigparm);
 507
 508    if(verbose>0){
 509      fprintf(dump, "coefs examined\n");
 510      tab_display(tp, dump);
 511      fflush(dump);
 512    }
 513
 514    if(pivi < nligne) goto pirouette;
 515    /* There is a row whose coefficients are negative */
 516    compa_test(tp, context, ni, nvar, nparm, nc);
 517    if(verbose>0){
 518      fprintf(dump, "compatibility tested\n");
 519      tab_display(tp, dump);
 520      fflush(dump);
 521    }
 522
 523    pivi = chercher(tp, Minus, nligne);
 524    if(pivi < nligne) goto pirouette;
 525    /* The compatibility test has found a negative row */
 526    pivi = chercher(tp, Critic, nligne);
 527    if(pivi >= nligne)pivi = chercher(tp, Unknown, nligne);
 528    /* Here, the problem tree splits        */
 529    if(pivi < nligne) {
 530      Tableau * ntp;
 531      Entier com_dem;
 532      struct high_water_mark q;
 533      if(nc >= context->height) {
 534        dcw = mpz_sizeinbase(context->determinant,2);
 535        dch = 2 * dcw + 1;
 536        context = expanser(context, 0, nc, nparm+1, 0, dch, dcw);
 537      }
 538      if(nparm >= MAXPARM) {
 539        fprintf(stdout, "Too much parameters : %d\n", nparm);
 540        exit(2);
 541      }
 542      q = tab_hwm();
 543      if(verbose>0)
 544        fprintf(stdout,"profondeur %d %lx\n", profondeur, q.top);
 545      ntp = expanser(tp, nvar, ni, ncol, 0, 0, 0);
 546      fflush(stdout);
 547      sol_if();
 548      sol_forme(nparm+1);
 549      mpz_init_set_ui(com_dem, 0);
 550      for(j = 0; j<nparm; j++) {
 551        mpz_set(discr[j], Index(tp, pivi, j + nvar +1));
 552        mpz_gcd(com_dem, com_dem, discr[j]);
 553      }
 554      mpz_set(discr[nparm], Index(tp, pivi, nvar));
 555      mpz_gcd(com_dem, com_dem, discr[nparm]);
 556      for(j = 0; j<=nparm; j++) {
 557        mpz_divexact(discr[j], discr[j], com_dem);
 558        mpz_set(Index(context, nc, j), discr[j]);
 559        sol_val(discr[j], UN);
 560      }
 561      mpz_clear(com_dem);
 562      Flag(context, nc) = Unknown;
 563      mpz_set(Denom(context, nc), UN);
 564      Flag(ntp, pivi) = Plus;
 565      profondeur++;
 566      fflush(stdout);
 567      if(verbose > 0) fflush(dump);
 568      traiter(ntp, context, iq, nvar, nparm, ni, nc+1, bigparm);
 569      profondeur--;
 570      tab_reset(q);
 571      if(verbose>0)
 572        fprintf(stdout, "descente %d %lx\n", profondeur, tab_hwm().top);
 573      for(j = 0; j<nparm; j++)
 574        mpz_neg(Index(context, nc, j), Index(context, nc, j));
 575      mpz_add_ui(Index(context, nc, nparm), Index(context, nc, nparm), 1);
 576      mpz_neg(Index(context, nc, nparm), Index(context, nc, nparm));
 577      Flag(tp, pivi) = Minus;
 578      mpz_set(Denom(context, nc), UN);
 579      nc++;
 580      goto pirouette;
 581    }
 582 /* Here, all rows are positive. Do we need an integral solution?      */
 583    if(!iq) {
 584      solution(tp, nvar, nparm);
 585      break;
 586    }
 587 /* Yes we do! */
 588    pivi = integrer(&tp, &context, &nvar, &nparm, &ni, &nc);
 589    if(pivi > 0) goto pirouette;
 590                     /* A cut has been inserted and is always negative */
 591 /* Here, either there is an integral solution, */
 592    if(pivi == 0) solution(tp, nvar, nparm);
 593 /* or no solution exists */
 594    else sol_nil();
 595    break;
 596
 597 /* Here, a negative row has been found. The call to <<pivoter>> executes
 598       a pivoting step                                                 */
 599
 600 pirouette :
 601      if(pivoter(tp, pivi, nvar, nparm, ni, iq) < 0) {
 602        sol_nil();
 603        break;
 604      }
 605  }
 606 /* Danger : a premature return would induce memory leaks   */
 607  tab_reset(x);
 608  for(i=0; i<MAXPARM; i++)
 609    mpz_clear(discr[i]);
 610  return;
 611 }
 612
 613
 614
 615
 616
 617
 618