/* ========================================================================= */ /* SIMPLEXE for integer variables */ /* ALL-INTEGER METHODS */ /* Jean Claude SOGNO */ /* Projet CHLOE -- INRIA ROCQUENCOURT */ /* Juin 1994 */ /* ========================================================================= */ /* ========================================================================= */ /* Duong NGUYEN QUE */ /* Adaption to abstract computation: janusvalue */ /* CRI-ENSMP */ /* ========================================================================= */ #ifdef HAVE_CONFIG_H #include "config.h" #endif #include "arithmetique.h" #include #include "assert.h" #include "iproblem.h" #include "iabrev.h" #include "rproblem.h" /* ============================================================== */ /* PROCEDURES for printing computation information */ /* ============================================================== */ static int xbug(Pproblem XX,int nu) { printf("DIRECT xBug=%d Pmeth=%d Pmet2=%d Pmet3=%d\\\\\n", nu,PMETH,PMET2,PMET3); if (VISU||VIS2) fprintf(FTRACE,"xBug=%d\\\\\n",nu); return(VRBUG); } static int xdeb(Pproblem XX,int nu) { if (VISU||VIS2) fprintf(FTRACE,"xDeb=%d\\\\\n",nu); return(VRDEB); } static void vzphase(Pproblem XX,int n) /*****************************/ { if (VISU>=ZVS) { fprintf(FTRACE,"{\\bf{*} PHASE: "); if (n==1) fprintf(FTRACE,"PRE-PROCESSING PHASE"); else if (n==2)fprintf(FTRACE,"DUAL SIMPLEX"); else if (n==3)fprintf(FTRACE,"PRIMAL SIMPLEX"); fprintf(FTRACE,"}\\\\\n"); } } static void vzstep(Pproblem XX,int n) /*****************************/ { if (VISU>=ZVS) { fprintf(FTRACE,"{\\bf STEP %d: ",n); if (n==1) fprintf(FTRACE,"EQUALITIES ELIMINATION"); else if (n==2)fprintf(FTRACE,"INEQUALITIES DIVISION"); else if (n==3)fprintf(FTRACE,"FOURIER-MOTZKIN"); else if (n==4)fprintf(FTRACE,"NON-NEGATIVE VARIABLES"); else if (n==5)fprintf(FTRACE,"FOURIER-MOTZKIN on CONSTRAINED VARIABLES"); else if (n==6) fprintf(FTRACE,"REDUNDANT INEQUALITIES"); fprintf(FTRACE,"}\\\\\n"); } } static void vznowstep(Pproblem XX) /*****************************/ { if (MAJSTEP) return; vzstep(XX,NSTEP); MAJSTEP=1; } static void vznewstep(Pproblem XX,int n) /*****************************/ { NSTEP=n; MAJSTEP=0; } void tableau_janus(Pinitpb II, Pproblem XX) /**************************/ { int i,j,i00; i00=1 ; if (NV>100) return; fprintf(FTRACE,"%%%d %d %d %d\n",NV,MC,NP,i00); /*****fprintf(FTRACE,"%%%d %d %d %d\n",XX->nvar,XX->MCONTR,XX->NVPOS,i00);**/ for (i=i00 ; i<=MC ; i++ ) { fprintf(FTRACE,"%%%4d",E(i)); for (j=1 ; j <= NV ; j++) fprint_Value(FTRACE,A(i,j)); fprint_Value(FTRACE,D(i)); fprintf(FTRACE," \n"); /*fprintf(FTRACE,"%%%4d",E(i)); for (j=1 ; j <= NV ; j++) fprintf(FTRACE," %2d",A(i,j)); fprintf(FTRACE," %3d\n",D(i));DN*/ } /*fprintf(FTRACE," \\\\\n");*/ } /*DN note: if use janus.c then not static DN*/ /* static void tableauk(Pproblem XX,int j1,int j2) */ /* { int i,j,i00; i00=1 ; */ /* if (j1>NX) {fprintf(FTRACE," --- pas de colonne %d ----\n",j1); return; }; */ /* if (j2>NX) {fprintf(FTRACE," --- pas de colonne %d ----\n",j2); return; }; */ /* fprintf(FTRACE,"%d %d %d %d\n",NX,MX,NX,i00); */ /* for (i=i00 ; i<=MX ; i++ ) */ /* { fprintf(FTRACE,"1 "); */ /* if (j2>0) */ /* { for (j=1 ; j <= NX ; j++) */ /* { if (j==j1) fprint_Value(FTRACE,AK(i,j2)); */ /* else if (j==j2) fprint_Value(FTRACE,AK(i,j1)); */ /* else fprint_Value(FTRACE,AK(i,j)); */ /* /\*if (j==j1) fprintf(FTRACE," %2d",AK(i,j2)); */ /* else if (j==j2) fprintf(FTRACE," %2d",AK(i,j1)); */ /* else fprintf(FTRACE," %2d",AK(i,j));DN*\/ */ /* } */ /* } */ /* else */ /* { for (j=1 ; j <= NX ; j++) if (j!=j1) fprint_Value(FTRACE,AK(i,j)); */ /* if (j1>0 && (j1 <= NX)) fprint_Value(FTRACE,AK(i,j1)); */ /* /\*for (j=1 ; j <= NX ; j++) if (j!=j1) fprintf(FTRACE," %2d",AK(i,j)); */ /* if (j1>0 && (j1 <= NX)) fprintf(FTRACE," %2d",AK(i,j1));DN*\/ */ /* } */ /* fprint_Value(FTRACE,DK(i)); fprintf(FTRACE," \n"); */ /* /\*fprintf(FTRACE," %3d\n",DK(i));DN*\/ */ /* } */ /* } */ static void symbol(Pproblem XX,int v) /*****************************/ { if (v>NUMAX) fprintf(FTRACE,"w") ; else if (v<0) fprintf(FTRACE,"v") ; else if (freevariable(XX,v)) fprintf(FTRACE,"x") ; else fprintf(FTRACE,"y") ; fprintf(FTRACE,"_{%d}",abs(v)); } static void symbold(Pproblem XX,int v) /****************************/ { fprintf(FTRACE,"$");symbol(XX,v);fprintf(FTRACE,"$"); } /* static void symboldn(Pproblem XX,int v) /\***************************\/ */ /* { symbold(XX,v);fprintf(FTRACE,"\\\\\n"); */ /* } */ static void vzcut(Pproblem XX,int i,int divisor,int i2) { int line,j ; if (NX>16) fprintf(FTRACE,"{\\scriptsize\n"); else if (NX>13) fprintf(FTRACE,"{\\footnotesize\n"); fprintf(FTRACE,"\\begin{center}\n"); fprintf(FTRACE,"\\fbox{$\n"); fprintf(FTRACE," \\begin{array}{rrrrrrrrrrrrrrrrrrrrrrrrr}\n"); for (line=1; line<=2; line++) { int s,nk;//unused ab Value k;//DN int fleche,carre,carre2; fleche=carre=carre2=0; s=0; nk=0; /*if (line==1) symbol(XX,G(i)); else fprintf(FTRACE," cut") ;*/ if (line==1) symbol(XX,G(i)); else symbol(XX,G(i2)); fprintf(FTRACE,":") ; for ( j=1 ; j <= NX ; j++) /*j0*/ { int boite; boite=0; if (value_assign(k,AK(i,j))) ++nk; if (nk==17) { nk=0; fprintf(FTRACE," \\\\ & \n") ; } if (( NX <= 24) || value_notzero_p(k) ) fprintf(FTRACE," &") ; if value_notzero_p(k) { if ((value_neg_p(k)) && (line==1)) fprintf(FTRACE,"-"); else if (s) fprintf(FTRACE,"+"); if (line==1) { if (value_notone_p(value_abs(k))) fprint_Value(FTRACE,k) ; } else {// fprintf(FTRACE,"\\lfloor%d/%d\\rfloor ",k,divisor) ; fprintf(FTRACE,"\\lfloor");fprint_Value(FTRACE,k);fprintf(FTRACE,"/%d\\rfloor ",divisor) ;//DN_? divisor?? } symbol(XX,B(j)); s=1; } if (boite) { fprintf(FTRACE,"$}"); boite=0; } } if ((i13) fprintf(FTRACE,"}\n"); } static void igvoir3(Pproblem XX,int ia,int ib,int carre, int ip, int jp, int carre2, int ip2, int jp2, int fleche, int ifl) { int i,j ; /*print latex arrays * ak et dk between lines ia and ib */ if (NX>16) fprintf(FTRACE,"{\\scriptsize\n"); else if (NX>13) fprintf(FTRACE,"{\\footnotesize\n"); fprintf(FTRACE,"\\begin{center}\n"); fprintf(FTRACE,"\\fbox{$\n"); fprintf(FTRACE," \\begin{array}{rrrrrrrrrrrrrrrrrrrrrrrrr}\n"); for ( i=ia ; i <= ib ; i++) { int s,nk; Value k;//DN, unused ab s=0; nk=0; if (fleche && i==ifl) fprintf(FTRACE,"\\Rightarrow ") ; if (i24) fprintf(FTRACE," .....") ; else for (j=1; j <= NX; j++) fprintf(FTRACE," &.....") ; /*j0*/ fprintf(FTRACE," \\\\\n") ; } } ; fprintf(FTRACE," \\end{array} \n") ; fprintf(FTRACE,"$}\n"); fprintf(FTRACE,"\\end{center}\n"); if (NX>14) fprintf(FTRACE,"}\n"); } /* fin igvoir3 */ static void igvoir(Pproblem XX,int ia,int ib) /*print latex arrays */ { igvoir3(XX,ia,ib,0,0,0,0,0,0,0,0); } /* fin igvoir */ static void w1voir(Pproblem XX,int ix) { igvoir(XX,ix,ix); /******* print line ak (ix) + dk (ix) */ } /*static void wvoir(Pproblem XX)*/ /**** display arrays ak and dk */ static void wvoir(Pproblem XX) /**** acces pour mise au point janus */ { if (NV>100) { int i,j; Value cab,cmax;//DN coeef max cmax = VALUE_ZERO; for ( i=1 ; i <= MX ; i++) for ( j=1 ; j <= NX ; j++) if (value_assign(cab,value_abs(AK(i,j)))) if value_gt(cab,cmax) value_assign(cmax,cab); {fprintf(FTRACE,"COEFFICIENT MAXIMUM=");fprint_Value(FTRACE,cmax);fprintf(FTRACE," \n");} return; } igvoir(XX,1,MX); /*i0*/ } /* static void stepvoir(Pproblem XX) /\**** acces pour mise au point janus *\/ */ /* { */ /* if (NV>100) { */ /* int i,j; */ /* Value cab,cmax;//DN coeef max */ /* cmax = VALUE_ZERO; */ /* for ( i=1 ; i <= MX ; i++) */ /* for ( j=1 ; j <= NX ; j++) */ /* if (value_assign(cab,value_abs(AK(i,j)))) */ /* if value_gt(cab,cmax) value_assign(cmax,cab); */ /* {fprintf(FTRACE,"COEFFICIENT MAXIMUM=");fprint_Value(FTRACE,cmax);fprintf(FTRACE," \n");} */ /* return; */ /* } */ /* wvoir(XX); /\*i0*\/ */ /* } */ static void vzemptycol(Pproblem XX,int j) { fprintf(FTRACE,"empty column, variable $"); symbol(XX,B(j)); fprintf(FTRACE,"$\\\\\n"); } static void vzredundant(Pproblem XX,int i,int i2) { if (VISU>=ZVS) vznowstep(XX); if(VISU>=ZVR1) { if(i2<0) fprintf(FTRACE," empty inequality $y_{%d}$\\\\",G(i)); else if(i2==0) fprintf(FTRACE," useless inequality $y_{%d}$\\\\",G(i)); else fprintf(FTRACE, " redundant inequality $y_{%d}$ compared with $y_{%d}$\\\\\n",G(i),G(i2)); if(VISU>=ZVR4) igvoir3(XX,1,MX,0,0,0,0,0,0,1,i); } } static void vzgcd(Pproblem XX,int i, int gcd, int viz)//DN int gcd or Value gcd??? { if (VISU>=viz) fprintf(FTRACE, "* Polyedron proved empty by GCD test on ligne%3d gcd=%2d\\\\\n",i,gcd); if (VISU>=viz+2) wvoir(XX); } static void vzextgcd(Pproblem XX,int a1,int a2,int u1,int u2,int u3, int zvx) { if (VISU=ZVU1+2) { if (VISU=ZVU1) { fprintf(FTRACE,"unimodular change $"); symbol(XX,old1); fprintf(FTRACE,", "); symbol(XX,old2); fprintf(FTRACE," \\Rightarrow "); symbol(XX,new1); fprintf(FTRACE,", "); symbol(XX,new2); fprintf(FTRACE,"$\\\\\n"); } if (VISU>=ZVU1+5) wvoir(XX); } /* ======================================================================= */ /* FOR REAL SIMPLEX */ /* ======================================================================= */ static int integrite(float vf) /* ***************************** */ { int vi ; float vfp, epss; epss=0.00001; /*epss=0.000001;*/ if (vf<0) vfp= -vf ; else vfp= vf ; vi=vfp+epss ; if (vfp-vi0) testdessusdessous(XX,-vf); } void testsdessusdessous(Pproblem XX) { testdessusdessous(XX,0); testdessusdessous(XX,7.6543); testdessusdessous(XX,4.000002); testdessusdessous(XX,4.000001); testdessusdessous(XX,1.999995); testdessusdessous(XX,1.999999); } in case of test - end *******************/ /************************* for real simplex ********************************/ static void integerrealline(Pproblem XX,struct rproblem *RR, int ii, int ir) { int j; RR->d[ir]=DK(ii); if (ir==0) RR->e[0]=XX->minimum ; else RR->e[ir]=1 ; for (j=1 ; j<=NX ; j++) RR->a[ir][j]=AK(ii,j); if (MSR>=7) RR->iname[NX+ir]=G(ii); } static void integertoreal(Pproblem XX,struct rproblem *RR,int mini) { int i,j; RR->nvar=NX; RR->mcontr=MX-IC1+1; RR->nvpos=RR->nvar; RR->copie=0; RR->ntrace=XX->ntrac3; RR->meth=0; if (PRECR) RR->testp=0; else RR->testp=3; RR->ftrace=FTRACE; RR->base=0; RR->cost=0; RR->rhs=0; for (i=IC1 ; i<=MX ; i++) integerrealline(XX,RR,i,i+1-IC1); if (ICOUT) integerrealline(XX,RR,ICOUT,0); else { RR->d[0]=0; if (mini==1) RR->e[0]=1 ; else {RR->e[0]= -1 ;fprintf(FTRACE,"creation real mini %d\\\\\n",mini); } for (j=1 ; j<=NX ; j++) RR->a[0][j]=0; } if (mini<0) RR->e[0]= -1 ; RR->state=0; if (MSR>=7) for (j=1 ; j<=NX ; j++) RR->iname[j]=B(j); } static void messrsimplex(Pproblem XX,struct rproblem *RR,int r, int ph) { if (!r) {fprintf(FTRACE," REEL-PHASE %d:%13.6f \n",ph,RR->x[0]); return; } if(r==VRVID)fprintf(FTRACE,"(EMPTY)"); else if(r==2)fprintf(FTRACE,"(INFINI)"); else if(r==3)fprintf(FTRACE,"INSUF"); else if(r==8)fprintf(FTRACE,"(EPSILON)"); fprintf(FTRACE,"ANOMALIE PHASE %d REEL %d it=%d\\\\\n",ph,r,RR->iter); } static float localnewrsimplex(Pproblem XX,struct rproblem *RR,int min) { int rsr,rd=0; float rrd =0,rrp; if (MSR==0) { integertoreal(XX,RR,min); RR->meth=1; rd=realsimplex(RR); rrd=RR->x[0]; /*fprintf(FTRACE,"Real dual simplex x[0]=%13.6f\\\\\n",vpara.x[0]);*/ } integertoreal(XX,RR,min); rsr=realsimplex(RR); if (rsr) messrsimplex(XX,RR,rsr,1); rrp=RR->x[0]; /*if (VW6) fprintf(FTRACE,"REEL-PRIMAL:%13.6f \n",rrp);*/ if (MSR==0) { if (rd!=rsr) fprintf(FTRACE,"DIVERGENCE DUAL:%d PRIMAL:%d\n",rd,rsr); else if (VW6) fprintf(FTRACE, "REEL DUAL:%13.6f PRIMAL:%13.6f ECART:%13.6f\n",rrd,rrp,rrd-rrp); } else if (VW6) fprintf(FTRACE,"REEL PRIMAL:%13.6f \n",rrp); /*XX->itemp=rsr; XX->ftemp=RR->x[0]; * revoir */ return(RR->x[0]); } static float localrsimplex(Pproblem XX,int min) { struct rproblem ARR; localnewrsimplex(XX,&ARR,min); return(ARR.x[0]); } static void copybound(Pproblem XX,Pproblem YY) { int i; for (i=1; i<= XX->numax; i++) value_assign(YY->ibound[i],XX->ibound[i]); for (i=1; i<= XX->numax; i++) value_assign(YY->ilbound[i],XX->ilbound[i]); } static void vid(Pproblem XX) { XX->state= -1; } static void vidr(struct rproblem *RR) { RR->state= -1; } /*static int copystructure(Pproblem XX,Pproblem YY)*/ /*int copystructure(Pproblem XX,Pproblem YY)*/ static int copystructure(Pproblem XX,Pproblem YY) { int i,j; if (XX==YY) return(1); if (XX->state<0) return(2); YY->state=XX->state; YY->nvar=XX->nvar; YY->mcontr=XX->mcontr; YY->nvpos=XX->nvpos; YY->ftrace=XX->ftrace; YY->dyn=XX->dyn; YY->ntrace=XX->ntrace; YY->ntrac2=XX->ntrac2; YY->ntrac3=XX->ntrac3; YY->fourier=XX->fourier; YY->varc=XX->varc; YY->forcer=XX->forcer; YY->meth=XX->meth; YY->met2=XX->met2; YY->met3=XX->met3; YY->met4=XX->met4; YY->met5=XX->met5; YY->met6=XX->met6; YY->met7=XX->met7; YY->met8=XX->met8; YY->critermax=XX->critermax; YY->remove=XX->remove; YY->turbo=XX->turbo; YY->choixpiv=XX->choixpiv; YY->choixprim=XX->choixprim; YY->negal=XX->negal; YY->icout=XX->icout; YY->minimum=XX->minimum; YY->mx=XX->mx; YY->nx=XX->nx; YY->tmax=XX->tmax; YY->numero=XX->numero; YY->numax=XX->numax; YY->lastfree=XX->lastfree; YY->niter=XX->niter; /*YY->result=XX->result;*/ /*YY->jinfini=XX->jinfini;*/ /*YY->nub=XX->nub;*/ YY->ic1=XX->ic1; YY->ntp=XX->ntp; YY->vdum=XX->vdum; for (i=1; i<= XX->mx; i++) { for (j=1; j<= XX->nx; j++) value_assign(YY->ak[i][j],XX->ak[i][j]); value_assign(YY->dk[i],XX->dk[i]); YY->g[i]=XX->g[i]; } for (j=1; j<= XX->nx; j++) value_assign(YY->b[j],XX->b[j]); if (XX->state) copybound(XX,YY); return(0); } /* ========================================================================= */ /* PROCEDURES FOR INFORMATION ABOUT VARIABLES AND SYSTEM */ /* ========================================================================= */ static int janus_egalite(Pproblem XX,int i) /*********** deplacer *****************/ { return(i>MX-NEGAL); } /* static int costvariable(Pproblem XX) /\**************************\/ */ /* { return(NV+1); */ /* } */ static int freevariable(Pproblem XX,int v) /**************************/ { return((v>NP && v<=NV) || v<0); } //DNstatic int freevariable(Pproblem XX,Value v) /**************************/ //{ return((value_gt(v,int_to_value(NP)) && value_le(v,int_to_value(NV)) || value_neg_p(v)); //} static int cutvariable(Pproblem XX, int v) /**************************/ { return(v>NUMAX); } static int fluidvariable(Pproblem XX, int v) /******* free variable */ { return( freevariable(XX,v) || cutvariable(XX,v) ); /* or cut variable */ } static int cutline(Pproblem XX,int i) /*******************************/ { return(cutvariable(XX,G(i))); /* is the constraint a cut inequality ? */ } static int fluidline(Pproblem XX,int i) /*****************************/ { return(fluidvariable(XX,G(i))); } /* static int cutcolumn(Pproblem XX,int j) /\*****************************\/ */ /* { return(cutvariable(XX,B(j))); */ /* } */ /* static int fluidcolumn(Pproblem XX, int j) /\**************************\/ */ /* { return(fluidvariable(XX,B(j))); */ /* } */ static int freecolumn(Pproblem XX,int j) /****************************/ { return(freevariable(XX,B(j))); } //DNstatic int emptylhs(Pproblem XX,int i) /*****************************/ //{ int j; for (j=1; j<=NX; j++) if (AK(i,j)) return(0); // return(1); //} static int emptylhs(Pproblem XX,int i) /*****************************/ { int j; for (j=1; j<=NX; j++) if value_notzero_p(AK(i,j)) return(0);//DN return(1); } static int useless(Pproblem XX,int i) /*****************************/ { int j,w; Value a; w=1; if value_neg_p(DK(i)) return(0);//DN for (j=1; j<=NX; j++) if (value_assign(a,AK(i,j))) { if value_pos_p(a>0) return(0); if (cutvariable(XX,B(j))) w=2; } return(w); } static int presence(Pproblem XX,int iv) { int i,j; XX->cu=0; XX->lu=0; for (i=IC1 ; i<=MX ; i++) if (G(i)==iv) {XX->lu=i; return(-i);} for (j=1 ; j<=NX ; j++) if (B(j)==iv) return(XX->cu=j); return(0); } static int newfreevariable(Pproblem XX) /**************************/ { return(--LASTFREE); } /* ======================================================================== */ /* PROCEDURES for data overflow test */ /* ======================================================================== */ //DNstatic int correctm(Pproblem XX,int k1,int k2) //{ int k ; k=k1*k2 ; /*examen si multiplication termes non nuls sans overflow*/ // if (k1>0 && k2>0 || k1<0 && k2<0) // { if (k>0) return(k); } // else if (k<0) return(k); // if (VISU>=ZVO) // { fprintf(FTRACE,"t:%4d k1=%12d,k2=%12d,k=%12d mult deb*\n",NITER,k1,k2,k); // } // return(0); //} //static int corrects(Pproblem XX,int k1,int k2) /*provisoire*/ //{ int k; k=k1+k2; /* examen si somme termes dont un non nul sans overflow */ // if (k1>=0 && k2>=0 && k<=0 || k1<=0 && k2<=0 && k>=0) // { if(VISU>=ZVO)fprintf(FTRACE,"t:%4d k1=%12d,k2=%12d,k=%12d somm deb***\n", // NITER,k1,k2,k); return(0); // } // return(1); //} static Value dn_multiply(Value v1, Value v2) { Value v; if (value_zero_p(v1) || value_zero_p(v2)) return(VALUE_ZERO); v = value_direct_multiply(v1,v2); if value_eq(v1,value_div(v,v2)) return(v); else { fprintf(stderr,"\nDNDNDN JANUS WARNING, multiplication overflow"); assert(false); return VALUE_NAN; } } static Value correctm(Pproblem XX,Value k1,Value k2) { Value k; if (value_notzero_p(k1) && value_notzero_p(k2)) {// must not be zero k = value_direct_multiply(k1,k2); if value_eq(k1,value_div(k,k2)) return(k); } if (VISU>=ZVO) { fprintf(FTRACE,"\nDNDNDN JANUS WARNING: correctm overflow"); assert(false); } return(0);//should check parameter notzero before use, so return zero means overflow } static int corrects(Pproblem XX,Value k1,Value k2) { Value k; k=value_plus(k1,k2); if value_eq(k1,value_minus(k,k2)) return(1); if (VISU>=ZVO) { assert(false); fprintf(FTRACE,"\nDNDNDN JANUS WARNING: corrects overflow"); } return(0);//means overflow } /* ======================================================================== */ /* PIVOTING PROCEDURE */ /* ======================================================================== */ static int ipivotage2(Pproblem XX,int i1,int j1,int ipiv1,int ipiv2, int vtu) /* La matrice ak et les seconds membres dk sont entiers Le pivot est egal a 1 ou -1. le pivotage de la matrice s'effectue a partir de la ligne i1 et de la colonne j1. Les numeros de la variable d'ecart (basique) de la ligne i1, et de la variable independante (hors-base) de la colonne j1, contenus respectivement dans les tableaux g et b, sont permutes. Dans le cas de l'elimination d'une egalite, l'operation equivaut a chasser de la base une variable d'ecart identiquement nulle */ { int i,j; /* variables travail */ Value pv,p2,p3,pd,pp,dn_tmp; i = 0;j=0; value_assign(pv,VALUE_ZERO);value_assign(p2,VALUE_ZERO);value_assign(p3,VALUE_ZERO); value_assign(pd,VALUE_ZERO);value_assign(pp,VALUE_ZERO);value_assign(dn_tmp,VALUE_ZERO); /*fprintf(stderr,"%d %d %d %d %d \n",i1,j1,ipiv1,ipiv2,vtu);*///DNDB if (NITER >=TMAX) return(VRESULT=VRINS); if (DYN) dynam(XX,0); NITER +=1; VRESULT=VRFIN; if (vtu) value_assign(pv,XX->tturb[i1][j1]); else { value_assign(pv,AK(i1,j1)); j=B(j1); B(j1)=G(i1); G(i1)=j ; }//DN /*fprintf(stderr,"%d %d ",i,j);fprint_Value(stderr,pv);fprintf(stderr," "); fprint_Value(stderr,p2);fprintf(stderr," ");fprint_Value(stderr,p3);fprintf(stderr," "); fprint_Value(stderr,pd);fprintf(stderr," ");fprint_Value(stderr,pp);fprintf(stderr," "); fprint_Value(stderr,dn_tmp);fprintf(stderr," \n");*///DNDB if (value_notone_p(value_abs(pv))) { /*if(VISU)*/ fprintf(FTRACE,"*** bug,pivot different +-1:"); fprint_Value(FTRACE,pv);fprintf(FTRACE,"\n"); return(VRESULT=VRBUG); } if (XX->turbo) { if (vtu) { for (i=ipiv1; i<=ipiv2; i++) if (value_assign(p2,AK(i,j1))) { //p2 <> 0 if value_neg_p(pv) value_oppose(p2); //DK(i)-=XX->dturb[i1]*p2; if value_notzero_p(XX->dturb[i1]) value_assign(dn_tmp,dn_multiply(XX->dturb[i1],p2)); else value_assign(dn_tmp,VALUE_ZERO); value_substract(DK(i),dn_tmp);//DN for (j=1; j<=NX; j++) if (j==j1) value_assign(AK(i,j),value_uminus(p2)); else if (value_assign(p3,XX->tturb[i1][j])) { //AK(i,j)-=p2*p3; if (value_assign(dn_tmp,dn_multiply(p2,p3))) value_substract(AK(i,j),dn_tmp);//DN } } return(VRESULT=0); } for (i=ipiv1; i<=ipiv2; i++) if (i != i1 && (value_assign(p2,AK(i,j1)))) { if (value_neg_p(pv)) value_oppose(p2); //DK(i)-=DK(i1)*p2 ; if (value_assign(dn_tmp,dn_multiply(DK(i1),p2))) value_substract(DK(i1),dn_tmp); for (j=1; j<=NX; j++) if (j==j1) value_assign(AK(i,j),value_uminus(p2)); else if (value_assign(p3,AK(i1,j))) { //AK(i,j)-=p2*p3; if (value_assign(dn_tmp,dn_multiply(p2,p3))) value_substract(AK(i,j),dn_tmp);//DN } } if (value_mone_p(pv)) { value_oppose(DK(i1)); for (j=1; j<=NX; j++) if (j!=j1) value_oppose(AK(i1,j)); } } else { for (i=ipiv1; i<=ipiv2; i++) if ((value_assign(p2,AK(i,j1))) && i != i1) { if (value_neg_p(pv)) value_oppose(p2); if (value_assign(pd,DK(i1))) { if (value_assign(pp,dn_multiply(pd,p2))) { if (corrects(XX,DK(i),pp)==0) return(VRESULT=VROVF); value_substract(DK(i),pp); } else return(VRESULT=VROVF); } for (j=1 ; j<=NX ; j++) /*j0*/ if (j==j1) value_assign(AK(i,j),value_uminus(p2)); else if (value_notzero_p(AK(i1,j))) { if (value_zero_p(value_assign(p3,correctm(XX,p2,AK(i1,j))))) return(VRESULT=VROVF); if (corrects(XX,AK(i,j),value_uminus(p3))==0) return(VRESULT=VROVF); value_substract(AK(i,j),p3); } } if (value_mone_p(pv)) { value_oppose(DK(i1)); for (j=1 ; j<=NX ; j++) /*j0*/ if (j!=j1) value_oppose(AK(i1,j)) ; } } if (VISU>=ZVP1) { fprintf(FTRACE,"{\\bf iteration} %d: ",NITER); fprintf(FTRACE,"pivot(%3d,%3d), ",i1,j1); fprintf(FTRACE," variables $"); if (janus_egalite(XX,i1)) fprintf(FTRACE,"e_{%d}",B(j1)); else symbol(XX,B(j1)); fprintf(FTRACE,"$ - $"); symbol(XX,G(i1)); fprintf(FTRACE,"$\\\\\n"); if (VISU>=ZVP2) wvoir(XX) ; } /*fprintf(stderr,"ketthucipivotage2 :i %dj %d pv",i,j);fprint_Value(stderr,pv);fprintf(stderr," p2"); fprint_Value(stderr,p2);fprintf(stderr," p3");fprint_Value(stderr,p3);fprintf(stderr," pd"); fprint_Value(stderr,pd);fprintf(stderr," pp");fprint_Value(stderr,pp);fprintf(stderr," dn_tmp"); fprint_Value(stderr,dn_tmp);fprintf(stderr," \n");*///DNDB return(VRESULT) ; } static int ipivotage(Pproblem XX,int i1,int j1) { return ipivotage2(XX,i1,j1,1,MX,0); } static int redundant(Pproblem XX,int i1,int i2) /**********************/ { int j; Value a,s; value_assign(s,value_minus(DK(i1),DK(i2))); for (j=1; j<=NX; j++) if (value_assign(a,value_minus(AK(i2,j),AK(i1,j)))) { if (cutvariable(XX,B(j))) return(0); if (value_notzero_p(s)) { if (value_pos_p(s)) { if (value_neg_p(a)) return(0); } else if (value_pos_p(a)) return(0); } else value_assign(s,a); } if (value_pos_p(s)) return(i1); else return(i2); } static int satisfait(Pproblem XX) /****** for a system including only */ { int i; /* inequalities, does a basic solution is obvious? */ for (i=MX; i>=IC1; i--) if (value_neg_p(DK(i))) return(0); return(1) ; } /* ========================================================================= */ /* PROCEDURE FOR REPLACING A CONSTRAINED VARIABLE BY A FREE VARIABLE */ /* ========================================================================= */ static int makecolumnfree(Pproblem XX,int jx) /************************/ { if (freecolumn(XX,jx)) return(0); if (cutvariable(XX,B(jx))) B(jx)= newfreevariable(XX); else { int jj; if (++MX >MAXMX) return(VRDEB); G(MX)= B(jx) ; value_assign(DK(MX),VALUE_ZERO); for (jj=1 ; jj<=NX ; jj++) value_assign(AK(MX,jj),VALUE_ZERO); value_assign(AK(MX,jx),VALUE_MONE); B(jx) = newfreevariable(XX); if (VISU>=ZVVF) { fprintf(FTRACE, "constrained feature of column %2d is specified by an inequality:\n",jx); w1voir(XX,MX); } } return(0); } /* ========================================================================= */ /* PROCEDURES FOR HANDLING MATRIX */ /* ========================================================================= */ static void celimination(Pproblem XX,int j1) /*************************/ { if (j1 != NX) /* column j1 is eliminated (and replaced by column NX) */ { int i; B(j1)=B(NX) ; for (i=1 ; i<=MX ; i++) value_assign(AK(i,j1),AK(i,NX)) ; /*i0*/ } ; --NX ; } static void emptycolelim(Pproblem XX,int j1) /*************************/ { if (j1 != NX) /* empty column j1 is eliminated and replaced by column NX */ { int i; B(j1)=B(NX) ; for (i=1 ; i<=MX ; i++) /*i0*/ if (value_notzero_p(AK(i,NX))) value_assign(AK(i,j1),AK(i,NX)); } --NX ; } static void permutec(Pproblem XX,int j1, int j2) /*********************/ { if (j1 != j2) /* columns j1 and j2 are permuted */ { int i,dn_tmp; Value vi; dn_tmp=B(j1); B(j1)=B(j2) ; B(j2)=dn_tmp; for (i=1 ; i<=MX ; i++) /*i0*/ { value_assign(vi,AK(i,j1)); value_assign(AK(i,j1),AK(i,j2)); value_assign(AK(i,j2),vi); } } } static void permutel(Pproblem XX,int i1, int i2) /*********************/ { if (i1 != i2) /* lines i1 and i2 are permuted */ { int j,v; Value vi; v=G(i1) ; G(i1)=G(i2) ; G(i2)=v ; value_assign(vi,DK(i1)); value_assign(DK(i1),DK(i2)); value_assign(DK(i2),vi); for (j=1 ; j<=NX ; j++) /*j0*/ { value_assign(vi,AK(i1,j)); value_assign(AK(i1,j),AK(i2,j)); value_assign(AK(i2,j),vi); } } } static void oppositeline(Pproblem XX,int io) /************/ { int j ; value_oppose(DK(io)); for (j=1; j<=NX; j++) value_oppose(AK(io,j)) ; /*j0*/ } static void retirer(Pproblem XX,int i1) /******************************/ { if (i1 != MX) /* line i1 is removed (and replaced by line MX) */ { int j ; value_assign(DK(i1),DK(MX)); G(i1)=G(MX); for (j=1; j<=NX; j++) value_assign(AK(i1,j),AK(MX,j)); /*j0*/ } --MX ; } /* static void badretirerineq(Pproblem XX,int i1) /\***********************\/ */ /* { if (NEGAL) permutel(XX,i1,MX-NEGAL); retirer(XX,MX-NEGAL); */ /* } */ static void retirerineq(Pproblem XX,int i1) /***********************/ { permutel(XX,i1,MX-NEGAL); retirer(XX,MX-NEGAL); } static void razcut(Pproblem XX) /******************************/ { int i; for (i=MX; i>=IC1; i--) if (cutline(XX,i)) retirerineq(XX,i); } /* ======================================================================= */ /* PROCEDURES FOR FOURIER-MOTZKIN */ /* ======================================================================= */ static int possible(Pproblem XX,int j) { int i; /** Fourier-Motzkin should not spoil functions ********/ for (i=1 ; i=ZVS) vznowstep(XX); if (VISU>=ZVF1) { fprintf(FTRACE,"fourier-Motzkin 0 : variable $x_{%d}$ - inequalities:", B(j)); /*vvv*/ for (i=MX ; i>=1 ; i--) if (value_notzero_p(AK(i,j))) fprintf(FTRACE," $y_{%d}$ ",G(i)) ; if (VISU=IC1 ; i--) if (value_notzero_p(AK(i,j))) retirer(XX,i) ; emptycolelim(XX,j) ; if (VISU>=ZVF3) wvoir(XX); } /* fin fourier0 */ static void ajouter(Pproblem XX,Value sk, int ip1, int i) /**************/ { int jk; Value vt,dn_tmp; if (value_one_p(sk)) { for (jk=NX ; jk>=1 ; jk--) /*j0*/ if ( value_assign(vt,AK(ip1,jk))) value_addto(AK(i,jk),vt); value_addto(DK(i),DK(ip1)); } else { for (jk=NX ; jk>=1 ; jk--) /* if (jk !=j) si j transmis */ /*j0*/ if ( value_assign(vt,AK(ip1,jk))) { if (value_notzero_p(sk)) value_assign(dn_tmp,dn_multiply(vt,sk)); else value_assign(dn_tmp,VALUE_ZERO); value_addto(AK(i,jk), dn_tmp); } if (value_notzero_p(sk)&&value_notzero_p(DK(ip1))) value_assign(dn_tmp,dn_multiply(DK(ip1),sk)); else value_assign(dn_tmp,VALUE_ZERO); value_addto(DK(i),dn_tmp); //DK(i)+= DK(ip1) *sk ; } } static void fourier1(Pproblem XX, int pn, int j, int ip1) /************/ /* integer fourier elimination criterion and: */ /* a single positive coefficient or a single negative coefficient */ { int i,bj; Value sk,vt; bj=B(j); if (VISU>=ZVS) vznowstep(XX); /*if (pn==1)*/ if (value_one_p(value_assign(vt,AK(ip1,j)))) { for (i=MX ; i>=IC1 ; i--) if (value_neg_p(value_assign(sk,AK(i,j)))) ajouter(XX,value_uminus(sk),ip1,i);/* ligne ip1 sur ligne i */ } else for (i=MX ; i>=IC1 ; i--) if (value_pos_p(value_assign(sk,AK(i,j)))) ajouter(XX,sk,ip1,i); retirer(XX,ip1) ; emptycolelim(XX,j) ; if (VISU) { if (VISU>=ZVF1) { fprintf(FTRACE,"fourier-Motzkin ");fprint_Value(FTRACE,vt);fprintf(FTRACE," variable $x_{%d}$\\\\\n",bj);} if (VISU>=ZVF3) wvoir(XX); } } /* fin fourier1 */ static void fourier22(Pproblem XX, int j) /****************************/ /* integer fourier elimination criterion and: */ /* 2 positive coefficients and 2 negative coefficients */ { int ip1=0,ip2=0,im1=0,im2=0,i,jk; Value sk,rp1,rp2=0,rm1,rm2=0,dp1,dp2,dm1,dm2,ap1,ap2,am1,am2; if (VISU>=ZVS) vznowstep(XX); rp1=0 ; rm1=0; for (i=MX ; i>=IC1 ; i--) if (value_assign(sk,AK(i,j))){ if (value_pos_p(sk)) { if (value_zero_p(rp1)) { value_assign(rp1,sk) ; ip1=i; } else { value_assign(rp2,sk) ; ip2=i; } } else { if (value_zero_p(rm1)) { value_assign(rm1,sk) ; im1=i; } else { value_assign(rm2,sk) ; im2=i; } }} if (VISU>=ZVF1) fprintf(FTRACE, /*vvv*/ "Fourier-Motzkin 2-2 variable $x_{%d}$ inequalities +: $y_{%d}$ $y_{%d}$ inequalities -: $y_{%d}$ $y_{%d}$\\\\\n",B(j),G(ip1),G(ip2),G(im1),G(im2)); value_assign(dp1,DK(ip1)); value_assign(dp2,DK(ip2)); value_assign(dm1,DK(im1)); value_assign(dm2,DK(im2)); value_assign(DK(im1),value_minus(dn_multiply(rp1,dm1),dn_multiply(rm1,dp1))); value_assign(DK(im2),value_minus(dn_multiply(rp2,dm1),dn_multiply(rm1,dp2))); value_assign(DK(ip1),value_minus(dn_multiply(rp1,dm2),dn_multiply(rm2,dp1))); value_assign(DK(ip2),value_minus(dn_multiply(rp2,dm2),dn_multiply(rm2,dp2))); //DK(im1) = rp1*dm1 - rm1*dp1; DK(im2) = rp2*dm1 - rm1*dp2; //DK(ip1) = rp1*dm2 - rm2*dp1; DK(ip2) = rp2*dm2 - rm2*dp2; for (jk=NX ; jk>=1 ; jk--) /*j0*/ { value_assign(ap1,AK(ip1,jk)); value_assign(ap2,AK(ip2,jk)); value_assign(am1,AK(im1,jk)); value_assign(am2,AK(im2,jk)); value_assign(AK(im1,jk),value_minus(dn_multiply(rp1,am1),dn_multiply(rm1,ap1))); value_assign(AK(im2,jk),value_minus(dn_multiply(rp2,am1),dn_multiply(rm1,ap2))); value_assign(AK(ip1,jk),value_minus(dn_multiply(rp1,am2),dn_multiply(rm2,ap1))); value_assign(AK(ip2,jk),value_minus(dn_multiply(rp2,am2),dn_multiply(rm2,ap2))); //AK(im1,jk),rp1*am1 - rm1*ap1; AK(im2,jk) = rp2*am1 - rm1*ap2; // AK(ip1,jk) = rp1*am2 - rm2*ap1; AK(ip2,jk) = rp2*am2 - rm2*ap2; } if (VISU>=ZVF2) wvoir(XX); emptycolelim(XX,j); if (VISU>=ZVF3) wvoir(XX); } /* fin fourier22 */ /* ======================================================================= */ /* PROCEDURES FOR UNIMODULAR CHANGES OF VARIABLES */ /* ======================================================================= */ //DNstatic int pgcd2(int z1, int z2) /********** GCD of two numbers *************/ //{ int p,v,r; p=abs(z1); v=abs(z2); // while (v!=0) // { r=p-(p/v)*v; p=v; v=r; // } /*fpprintf(FTRACE,"pgcd2(%3d,%3d)=%3d:\n",z1,z2,p);*/ // return(p); /* retour 0 si z1==z2==0 */ //} static Value pgcd2(Value z1, Value z2) /********** GCD of two numbers *************/ { Value p,v,r; value_assign(p,value_abs(z1)); value_assign(v,value_abs(z2)); while (value_notzero_p(v)) { value_assign(r,value_minus(p,dn_multiply(value_div(p,v),v))); value_assign(p,v); value_assign(v,r); } return(p); /* retour 0 si z1==z2==0 */ } /*static int extgcd(Pproblem XX,int a1,int a2,int *pu1,int *pu2,int *pu3, int zvx) { int v1,v2,v3,t1,t2,t3,q; if (a1==0 || a2==0) return(1); *pu1=1 ; *pu2=0 ; *pu3=abs(a1) ; v1=0 ; v2=1 ; v3=abs(a2) ; while (v3!=0) { q=*pu3/v3 ; t1= *pu1 -v1*q ; t2= *pu2-v2*q ; t3=*pu3-v3*q ; *pu1=v1 ; *pu2=v2 ; *pu3=v3 ; v1=t1 ; v2=t2 ; v3=t3 ; } if (a1<0) *pu1= - *pu1 ; if (a2<0) *pu2= - *pu2 ; if(VISU>=zvx) vzextgcd(XX, a1, a2, *pu1, *pu2, *pu3,zvx); return(0); }DN*/ static Value extgcd(Pproblem XX, Value a1, Value a2, Value *pu1,Value *pu2, Value *pu3, int zvx) { Value v1,v2,v3,t1,t2,t3,q; if (value_zero_p(a1) || value_zero_p(a2)) return(VALUE_ONE); value_assign(*pu1,VALUE_ONE); value_assign(*pu2,VALUE_ZERO); value_assign(*pu3,value_abs(a1)); value_assign(v1,VALUE_ZERO); value_assign(v2,VALUE_ONE); value_assign(v3,value_abs(a2)) ; while (value_notzero_p(v3)) { value_assign(q,value_div(*pu3,v3)) ; value_assign(t1,value_minus(*pu1,dn_multiply(v1,q))); value_assign(t2,value_minus(*pu2,dn_multiply(v2,q))); value_assign(t3,value_minus(*pu3,dn_multiply(v3,q))); value_assign(*pu1,v1); value_assign(*pu2,v2); value_assign(*pu3,v3); value_assign(v1,t1); value_assign(v2,t2); value_assign(v3,t3) ; } if (value_neg_p(a1)) value_oppose(*pu1) ; if (value_neg_p(a2)) value_oppose(*pu2) ; if(VISU>=zvx) vzextgcd(XX, a1, a2, *pu1, *pu2, *pu3,zvx); return(VALUE_ZERO); } static int jnsegment(Pproblem XX,int ik,int j1,int j2) { /*indique indice jk sur segment ligne ik tel que abs(coefficient)non nul, sinon renvoie -1 */ int j ; if (j1<=j2) { for (j=j1 ; j<=j2 ; j++) if (value_abs(AK(ik,j))) return(j); }else { for (j=j1 ; j>=j2 ; j--) if (value_abs(AK(ik,j))) return(j);} return(-1); } static int jsegment(Pproblem XX,int ik,int j1,int j2,Value p) { /* indique indice jk sur la ligne ik tel que abs(coefficient)==p , sinon renvoie -1 */ int j ; for (j=j1 ; j<=j2 ; j++) if (value_eq(value_abs(AK(ik,j)),p)) return(j); return(-1); } static Value pgcdsegment(Pproblem XX,int ik, int j1, int j2) /**********/ { int j; Value aa,p = VALUE_ZERO; for (j=j1 ; j<=j2 ; j++) if (value_assign(aa,AK(ik,j))) if (value_one_p(value_assign(p,pgcd2(p,aa)))) break; return(p); /* return: 0 if all coefficients null - otherwise gcd */ } static int changing2(Pproblem XX, int i1, int ja, int jb) /***********/ { int iz,old1,old2,new1,new2; Value za,zb, u1,u2,u3; value_assign(za,AK(i1,ja)) ; value_assign(zb,AK(i1,jb)); if (ja==jb) XBUG(7); old1=B(ja); old2=B(jb); new1=newfreevariable(XX); new2=newfreevariable(XX); if (VISU>=ZVU1+2) igvoir3(XX,1,MX,1,i1,ja,1,i1,jb,1,i1); else if (VISU>=ZVU1+1) igvoir3(XX,i1,i1,1,i1,ja,1,i1,jb,1,i1); if (extgcd(XX,za,zb,&u1,&u2,&u3,ZVU1+3)) XBUG(10); /*u3!=1 t97*/ // fprintf(stderr,"\n(in changing2: ");//DNDB for (iz=MX ; iz>=1 ; iz--) { Value aka,akb; value_assign(aka,AK(iz,ja)); value_assign(akb,AK(iz,jb)); // AK(iz,ja)= u1*aka + u2*akb ; //AK(iz,jb)= (za*akb - zb*aka)/u3 ; value_assign(AK(iz,ja),value_plus(dn_multiply(u1,aka),dn_multiply(u2,akb))) ; value_assign(AK(iz,jb),value_div(value_minus(dn_multiply(za,akb),dn_multiply(zb,aka)),u3)) ; // fprint_Value(stderr,AK(iz,ja));fprintf(stderr,"AKizja");//DNDB // fprint_Value(stderr,AK(iz,jb));fprintf(stderr,"AKizjb");//DNDB } //fprintf(stderr,"in changing2)\n");//DNDB B(ja)=new1; B(jb)=new2; if (VISU) vzunimod(XX,u1,u2,u3,za,zb, old1,old2,new1,new2); return(0); /* si pas d'incident */ } static int build0(Pproblem XX, int i1,int jj,int *pj1) /***************/ { int jn1=0,jn2=0 ; /* A unitary coefficient is built on equation i1 */ if (VISU>=ZVEV+1) fprintf(FTRACE, "Equation $e_{%d}$ does not include a unitary coefficient.\\\\\n",G(i1)); if ((jn1=jnsegment(XX,i1,NX,jj))<0) XBUG(5); /* jj,NX */ while (value_notone_p(AK(i1,jn1))) { jn2 = jnsegment(XX,i1,jn1-1,jj); if (jn2<0) XBUG(6); /*jn1+1,NX*/ if (makecolumnfree(XX,jn1)) XDEB(1); if (makecolumnfree(XX,jn2)) XDEB(2); if ((VRESULT=changing2(XX,i1,jn1,jn2))) return(VRESULT); } if (VISU>=ZVEV+1) fprintf(FTRACE,"Now, a unitary coefficient is available for a gaussian elimination:\n"); *pj1=jn1; return(0); } static int build1(Pproblem XX,int i1, int jj) /***********************/ { int vj=0,jn1=0,jn2=0;/*unitary or single coefficient built on inequality i1*/ Value delta= VALUE_ZERO; jn1=jnsegment(XX,i1,jj,NX);if (jn1<0) return(jn1);/*no non-zero coefficient*/ permutec(XX,jj,jn1); value_assign(delta,pgcdsegment(XX,i1,jj,NX)); if (VISU>=ZVNG) {fprintf(FTRACE,"Delta%3d=",i1);fprint_Value(FTRACE,delta);fprintf(FTRACE,"\\\\\n");} vj = jsegment(XX,i1,jj,NX,delta); /* recherche coefficient=delta */ if (vj>=0) /*on sait que valeur absolue d'un coefficient=pgcd*/ permutec(XX,jj,vj); else /* t98 2eme exemple RR1828 7.3.2 */ while (value_ne(value_abs(AK(i1,jj)),delta)) { jn2=jnsegment(XX,i1,jj+1,NX); if ((VRESULT=changing2(XX,i1,jj,jn2))) return(VRESULT); } return(0); } /* ======================================================================= */ /* PROCEDURES FOR CUT OPERATIONS */ /* ======================================================================= */ //DNstatic int cfloor(int v1, int v2) /********** v2 doit etre positif *********/ //{ if (v1==0) return(0); if (v1>0) return(v1/v2); return(((v1+1)/v2)-1); //} /* division reelle v1/v2=cfloor+reste * cfloor:entier * 0<=reste<1 */ static Value cfloor(Value v1, Value v2) /********** v2 doit etre positif *********/ { Value v; //fprintf(stderr,"(v1=");fprint_Value(stderr,v1);fprintf(stderr," ");//DNDB //fprintf(stderr,"v2=");fprint_Value(stderr,v2);fprintf(stderr," ");//DNDB if (value_zero_p(v1)) {/*fprintf(stderr,"cfloor VALUE_ZERO)");*/return(VALUE_ZERO); }//DNDB if (value_pos_p(v1)) {/* fprintf(stderr,"cfloor_pos ");fprint_Value(stderr,value_div(v1,v2));fprintf(stderr,")");*///DNDB return(value_div(v1,v2));} value_assign(v,value_div(value_plus(v1,VALUE_ONE),v2)); value_decrement(v); /*fprintf(stderr,"cfloor_neg ");fprint_Value(stderr,v);fprintf(stderr,")");*///DNDB return(v); } static int newcut(Pproblem XX) /*********/ { if (++MX >MAXMX) return(VRDEB); G(MX)= ++NUMERO ; return(0); } static void coupe(Pproblem XX, int i1, int i2, Value abpivot) /**********/ {int j2; /* fprintf(stderr,"\n(in coupe: abpivot = ");fprint_Value(stderr,abpivot);fprintf(stderr,"");*///DNDB value_assign(DK(i2),cfloor(DK(i1),abpivot)); /* fprintf(stderr," DKi2 =");fprint_Value(stderr,DK(i2));fprintf(stderr,"\n");*///DNDB for (j2=1 ; j2<=NX ; j2++) /*j0*/ { value_assign(AK(i2,j2),cfloor(AK(i1,j2),abpivot)); /* fprintf(stderr," AKi2j2 = ");fprint_Value(stderr,AK(i2,j2));fprintf(stderr,"\n");*///DNDB } // fprintf(stderr,"\nin coupe)\n");//DNDB if (VISU>=ZVC1) /*3*/ { if (VISU>=ZVC1+1) fprintf(FTRACE,"Let's build at "); /*zvc2 4*/ fprintf(FTRACE,"cut $"); symbol(XX,G(i2)); if(i1==i2) {fprintf(FTRACE,"$ (by "); fprint_Value(FTRACE,abpivot); fprintf(FTRACE,") from surrogate inequality");} else {fprintf(FTRACE,"$ (by "); fprint_Value(FTRACE,abpivot); fprintf(FTRACE,") from source row $y_{%d}$",G(i1));} if (VISU>=ZVC1+3) w1voir(XX,i2); /*zvc4 6*/ if (VISU>=ZVC1+2) /*zvc3 5*/ { fprintf(FTRACE," as follows:"); vzcut(XX,i1,abpivot,i2); } else fprintf(FTRACE,"\\\\\n"); /*suppression annulee 13juin 00 */ } } static int cutiteration(Pproblem XX, int i1, int j1, int turb) /******** coupe eventuelle + pivotage + retrait ************/ { int r=0,ix=0; Value abpivot =VALUE_ZERO; if (value_notone_p(value_assign(abpivot,value_abs(AK(i1,j1))))) { if (value_zero_p(abpivot)) { NUB=13; return(VRBUG); } if (cutline(XX,i1)) /* surrogate cut*/ { coupe(XX,i1,i1,abpivot); ix=i1; } else { if ((r=newcut(XX))) return(r); coupe(XX,i1,MX,abpivot); ix=MX; if (PMETH==9) /* PPP9 */ {localrsimplex(XX,1);fprintf(FTRACE,"(apres cut,avant piv)\\\\\n"); } } if (VISU>=ZVCP1) igvoir3(XX,1,MX,1,i1,j1,1,ix,j1,1,i1); } else { ix=i1; if (VISU>=ZVCP1) igvoir3(XX,1,MX,1,i1,j1,0,0,0,0,0); } if (turb) { int tu,j; tu= ++XX->nturb; XX->jturb[tu]=j1; value_assign(XX->dturb[tu],DK(ix)); if (XX->nturb>29) {fprintf(FTRACE,"turbo insuffisant %d\\\\\n",XX->nturb);return(VROVF);} for (j=1; j<=NX; j++) value_assign(XX->tturb[tu][j],AK(ix,j)); if ((r=ipivotage2(XX,ix,j1,IC1,MX,0))) return(r); } else if ((r=ipivotage(XX,ix,j1))) return(r); if (fluidline(XX,ix)) { if (VISU>=ZVCP2) wvoir(XX); retirer(XX,ix); } if (VISU>=ZVCP3) wvoir(XX); return(0); } /* ======================================================================= */ /* INEQUALITY DIVISION */ /* ======================================================================= */ //DNstatic int pgcd2in(int z1, int z2) /********** GCD of two numbers ************/ //{ int p,v,r; p=abs(z1); v=abs(z2); // while (v!=0) // { r=p-(p/v)*v; p=v; v=r; // } /*fpprintf(FTRACE,"pgcd2(%3d,%3d)=%3d:\n",z1,z2,p);*/ // return(p); /* retour 0 si z1==z2==0 */ //} static Value pgcd2in(Value z1, Value z2) /********** GCD of two numbers ************/ { Value p,v,r; value_assign(p,value_abs(z1)); value_assign(v,value_abs(z2)); while (value_notzero_p(v)) { value_assign(r,value_minus(p,dn_multiply(value_div(p,v),v))); value_assign(p,v);value_assign(v,r); } return(p); /* retour 0 si z1==z2==0 */ } static Value pgcdseg(Pproblem XX,int ik, int j1, int j2) /**********/ { int j; Value aa, p; value_assign(p,VALUE_ZERO); for (j=j1 ; j<=j2 ; j++) if (value_assign(aa,AK(ik,j))) if (value_one_p(value_assign(p,pgcd2in(p,aa)))) break; return(p); /* return: 0 if all coefficients null - otherwise gcd */ } static Value inequaldivision(Pproblem XX,int i) { int j; Value gcd; /* return 0: all coefficients are null * otherwise gcd */ if (value_zero_p(value_assign(gcd,pgcdseg(XX,i,1,NX)))) return(VALUE_ZERO); /* is LHS empty? */ else if (value_gt(gcd,VALUE_ONE)) { Value dp ; if (VISU>=ZVS) vznowstep(XX); if (VISU>=ZVI1) { if (VISU>=4) wvoir(XX); /* 2000 */ fprintf(FTRACE," line %d: inequality $",i); symbol(XX,G(i)); {fprintf(FTRACE," $ divided by gcd=");fprint_Value(FTRACE,gcd);fprintf(FTRACE,"\n");} if (VISU=1 ; j--) value_division(AK(i,j),gcd); /* division by GCD */ if (VISU>=ZVI3) w1voir(XX,i); } return(gcd); } /* ======================================================================= */ /* GCD TEST */ /* ======================================================================= */ static int gcdtest(Pproblem XX,int i, int *pj1, Value *ppivot, int viz) { int j; Value gcd; /* GCD test + division by GCD ** 0: fails * -1: empty polyedron */ if (value_zero_p(value_assign(gcd,pgcdsegment(XX,i,1,NX))))/* is LHS empty? */ { if (value_notzero_p(DK(i))) { if (VISU) vzgcd(XX,i,0,viz); return(-1); /*1*/ } } else if (value_gt(gcd,VALUE_ONE)) { Value dp ; if (VISU>=viz+1) { fprintf(FTRACE,"* ligne%3d gcd=",i);fprint_Value(FTRACE,gcd);fprintf(FTRACE,"\\\\\n");} /*2*/ if (value_notzero_p(DK(i))) /* GCD test */ { value_assign(dp,value_div(DK(i),gcd)); if (value_ne(dn_multiply(dp,gcd),DK(i))) { if (VISU) vzgcd(XX,i,gcd,viz); return(-1);} value_assign(DK(i),dp); } for (j=NX; j>=1; j--) value_division(AK(i,j),gcd); /* division by GCD */ for (j=NX; j>=1; j--) if (value_one_p(value_abs(AK(i,j)))) { *pj1=j ; value_assign(*ppivot,VALUE_ONE); break; } if (VISU>=viz+1){ fprintf(FTRACE,"Line %d was divided by gcd = ",i);fprint_Value(FTRACE,gcd);fprintf(FTRACE,"\\\\\n");} } return(0); } /* ======================================================================= */ /* BASIC ALL-INTEGER SIMPLEX PROCEDURES */ /* ======================================================================= */ static int dualentier(Pproblem XX) /********************************/ { int i,j,i1=0,j1,boul2,nn,nvn; Value teta,tet,pivot=VALUE_ZERO; if (VISU>=ZVS) { fprintf(FTRACE,"{*} STEP DUAL ALL.I. ,NX=%3d,MX=%3d,t=%4d,CHOIXPIV=%d ", NX,MX,NITER,CHOIXPIV); if (PMETH==MTDI2) fprintf(FTRACE," surrogate\\\\\n"); else fprintf(FTRACE," simple\\\\\n"); } XX->nturb=0; if (satisfait(XX)) return(VRFIN) ; for (;;) { int vg,r; vg=10000; value_assign(tet,VALUE_ZERO); nn=0 ; for (i=IC1 ; i<=MX ; i++) if (value_neg_p(value_assign(teta,DK(i)))) { boul2=1; ++nn ; for (j=1 ; j<=NX ; j++) if (value_neg_p(AK(i,j))) { boul2=0; break; } if (boul2) { if (VISU>=ZVDEND) wvoir(XX); return(VRVID); } if (CHOIXPIV==4) { if (vg>G(i)) { i1=i; value_assign(tet,teta); vg=G(i);} } else if (value_lt(teta,tet)) { i1=i; value_assign(tet,teta); } } if (value_zero_p(tet)) { if (XX->turbo>=3) if (XX->nturb) { int tu; for (tu=1 ; tu<=XX->nturb ; tu++) if ((r=ipivotage2(XX,tu,XX->jturb[tu],1,IC1-1,1))) return(r); XX->nturb=0; } if (VISU>=ZVDEND) wvoir(XX); return(VRFIN) ; /* available solution */ } nvn=0; for (j=1 ; j<=NX ; j++) if (value_neg_p(AK(i1,j))) { if (++nvn >1) break; j1=j ; value_assign(pivot, AK(i1,j)) ; } if (nvn >1) /* more than one possible pivot */ { if (CHOIXPIV==2||CHOIXPIV==4) { int comptj,compp; compp=0; for (j=1 ; j<=NX ; j++) if (AK(i1,j)<0) { comptj=0; for (i=IC1; i<=MX; i++) if ( DK(i) < 0) if (AK(i,j)<0) ++comptj; if (comptj>compp) { compp=comptj; j1=j; } } } else if (CHOIXPIV==3) { int comptj,compp; compp=0; for (j=1 ; j<=NX ; j++) if (value_neg_p(AK(i1,j))) { comptj=0; for (i=IC1 ; i<=MX ; i++) if ((value_neg_p(DK(i))) && (value_neg_p(AK(i,j)))) ++comptj; /*> > */ if (comptj>compp) { compp=comptj; j1=j; } } } else if (CHOIXPIV) { for (j=1 ; j<=NX ; j++) if (value_neg_p(AK(i1,j))) if (value_lt(pivot,AK(i1,j))) value_assign(pivot, AK(i1,j1=j)); } else for (j=1; j<=NX; j++) if (value_gt(pivot,AK(i1,j))) value_assign(pivot,AK(i1,j1=j)); } if (nn>1 && PMETH==MTDI2) /************ construction surrogate */ { int mx1; if (VISU>=ZVDS) fprintf(FTRACE, "surrogate inequality from %d inequalities\n",nn); mx1=MX ; if ((r=newcut(XX))) return(r); /*i1=MX;*/ value_assign(DK(MX),VALUE_ZERO); for (j=1 ; j<=NX ; j++) value_assign(AK(MX,j),VALUE_ZERO); for (i=IC1 ; i<=mx1 ; i++) if (value_neg_p(DK(i))) { value_addto(DK(MX),DK(i)) ; for (j=1 ; j<=NX ; j++) value_addto(AK(MX,j),AK(i,j)); } if (value_posz_p(DK(MX))) return(VROVF); if (VISU>=ZVDS) igvoir(XX,MX,MX); boul2=1; value_assign(pivot,VALUE_ZERO); /*** choix pivot surrogate ***/ if (CHOIXPIV) { for (j=1; j<=NX; j++) if (value_neg_p(AK(MX,j))) if (boul2 || value_lt(pivot,AK(MX,j))) { boul2=0 ; j1=j ; value_assign(pivot,AK(MX,j)) ; } } else for (j=1 ; j<=NX ; j++) if (value_neg_p(AK(MX,j))) if (boul2 || value_gt(pivot,AK(MX,j))) { boul2=0 ; j1=j ; value_assign(pivot, AK(MX,j)) ; } if (boul2) return(VRVID); i1=MX; } r=0; if (DYN) { if (VISU) { fprintf(FTRACE,"Dual, constraints:"); wvoir(XX);} r=dynam(XX,1); if (VISU) { if (r)fprintf(FTRACE,"sortie dynam Dual constraints: OUI\\\\\n"); else fprintf(FTRACE,"sortie dynam Dual constraints: non\\\\\n"); } } /*if (VISU) fprintf(FTRACE,"Etat sortie dyn=%d r=%d\\\\\n",DYN,r);*/ if (DYN && r) { if (VISU) fprintf(FTRACE,"unimodular change in DUAL\\\\\n"); if ((r=build0(XX,i1,1,&j1))) return(r); } if (XX->turbo>=2) { if ((r=cutiteration(XX,i1,j1,1))) return(r); if (XX->turbo==2) { int tu; tu=XX->nturb; if (tu!=1) {fprintf(FTRACE,"ARRET D URGENCE\\\\\n"); XBUG(37);} if ((r=ipivotage2(XX,tu,XX->jturb[tu],1,IC1-1,1))) return(r); XX->nturb=0; } } else if ((r=cutiteration(XX,i1,j1,0))) return(r); } } /* fin dualentier */ ; /* ======================================================================= */ static void razfreq(Pproblem XX) /*****/ { int i; for (i=1 ; i<=NUMAX ; i++) XX->frequency[i]=0; } /* static void voirfreq(Pproblem XX) /\*****\/ */ /* { int i; fprintf(FTRACE,"FREQUENCY:\\\\\n"); */ /* for (i=1 ; i<=NUMAX ; i++) if (XX->frequency[i]) */ /* fprintf(FTRACE,"frequence %d: %d\\\\\n",i,XX->frequency[i]); */ /* } */ static int iprimal(Pproblem XX,int minimum,int i2,int stopcout) /************ RUDIMENTARY PRIMAL ALL-INTEGER * cost line i2 ***********/ { int r,i,j,i1=0,j1=0,i3=0,bool1,lastit; float c1,c2,alpha=0,zeta ; Value pivot,lastdk,tet1,tet2; if (CRITERMAX) TMAX=(NX * CRITERMAX)+NITER; if (i2>=IC1 || !satisfait(XX)) { NUB=12; return(VRBUG);} if (VW6) fprintf(FTRACE,"Choix primal: %d (dual %d) iter %d $->$ %d", CHOIXPRIM,CHOIXPIV,NITER,TMAX); if (VISU) { vzphase(XX,3); if (minimum) fprintf(FTRACE," MIN"); else fprintf(FTRACE," MAX"); {fprintf(FTRACE,"IMIZATION function line %d (current value: ",i2); fprint_Value(FTRACE,value_uminus(DK(i2)));fprintf(FTRACE,")\\\\\n");} if (VISU>=ZVBR3) wvoir(XX) ; /* encore provisoire */ } if (NTPNTP ; j--) if (value_notzero_p(AK(i2,j))) { for (i=IC1 ; i<= MX ; i++) if (value_notzero_p(AK(i,j))) { NUB=14; return(VRBUG);} XX->jinfini=j; return(VRINF); } } value_assign(lastdk,DK(i2)); lastit=NITER; razfreq(XX); for (;;) { bool1=1; c1=0 ; value_assign(tet2,VALUE_ZERO); for (j=NX ; j>=1 ; j--) { if (value_neg_p(value_assign(tet1,(minimum) ? AK(i2,j) : value_uminus(AK(i2,j))))) { bool1=1; for (i=IC1 ; i<= MX ; i++) if (value_pos_p(AK(i,j))) { int selection; zeta = VALUE_TO_FLOAT(value_div(DK(i), AK(i,j))) ; if (CHOIXPRIM<=1) selection=(bool1||zeta=1.0)&&(zetaalpha)) || ((zeta<1.0)&&(alpha>=1.0)) ); if (selection) { bool1=0; alpha=zeta ; i3=i ; value_assign(pivot, AK(i,j)); } } if (bool1) { XX->jinfini=j; return(VRINF);} /*solution infinie*/ if (CHOIXPRIM==1) if (alpha<1.0) alpha=0; /* 10/2/98 */ if (((c2=alpha*(VALUE_TO_FLOAT(tet1)))frequency[G(i1)]; if ((r=cutiteration(XX,i1,j1,0))) return(r); /* PPP9 */ if (stopcout==9) return(stopcout); if (PMETH==9){localrsimplex(XX,1);fprintf(FTRACE,"(apres tou)\\\\\n");} if (value_ne(DK(i2),lastdk)) { if (VISU>=ZVPRI) {fprintf(FTRACE,"{\\bf iteration} %d improved cost value: ",NITER); fprint_Value(FTRACE,value_uminus(DK(i2)));fprintf(FTRACE,"\\\\\n");} value_assign(lastdk,DK(i2)); lastit=NITER; razfreq(XX); } if (stopcout && stopcout!=9) { if (value_posz_p(DK(i2))) return(10); for (i=i2+1; i< IC1 ; i++) if (value_posz_p(DK(i))) { razfreq(XX); return(11); } } } } /* fin iprimal ================================================ */ static void anonymouscopyline(Pproblem XX,int is, int id) /************/ { int j ; if (is==id) return; /* line is (source) is copied on line id */ value_assign(DK(id), DK(is)) ; for (j=1 ; j<=NX ; j++) value_assign(AK(id,j), AK(is,j)); /*j0*/ } static void copyline(Pproblem XX,int is, int id) /********************/ { if (is==id) return; /* line is (source) is copied on line id) */ anonymouscopyline(XX,is,id) ; G(id)=G(is) ; } /* ========================= PROCEDURES FOR FAST =========================== */ static void inhibit(Pproblem XX, int in) /*********/ { permutel(XX,in,IC1++); } static void remettre(Pproblem XX, int in) /*********/ { permutel(XX,in,IC1-1); IC1--; } static int remise(Pproblem XX, int icf) /*********/ { IC1=icf; retirer(XX,IC1); } static int fast(Pproblem XX,Pproblem ZZ) /*******************/ { int icf,iopt,r,i,j,exicou; exicou=ICOUT; icf=IC1 ; newcut(XX) ; inhibit(XX,MX); G(icf)=0; for (i=IC1; i<=MX; i++) if (value_neg_p(DK(i))) inhibit(XX,i); while (icfjinfini); wvoir(XX); if (value_mone_p(AK(icf,XX->jinfini))) { if (!iopt) XBUG(24); remettre(XX,iopt); ik=IC1; } else ik=icf; if ((r=cutiteration(XX,ik,XX->jinfini,0))) return(r); r=10; } ICOUT=exicou; if (value_neg_p(DK(icf)) && r!=11) /* no solution */ { remise(XX,icf); if (VISU>=ZVSAT3) { fprintf(FTRACE," fast vide\n"); wvoir(XX); } return(VRVID); } for (i=IC1-1 ; i>icf ; i--) if (value_posz_p(DK(i))) remettre(XX,i); if (VISU>=ZVSAT4) {fprintf(FTRACE,"STEP FAST IC1==%d\n",IC1); wvoir(XX);} } remise(XX,icf); if (VW2 || VISU>=ZVSAT4) { fprintf(FTRACE,"Feasibility proof (primal), cost = "); fprint_Value(FTRACE,value_uminus(DK(1)));fprintf(FTRACE,"\n"); vzlast(XX); wvoir(XX); } if (!satisfait(XX)) { NUB=27; return(VRBUG); } return(VRFIN); } /* ========================== fin fast =================================== */ static void anonymctol(Pproblem XX,int js, int id) /**********/ { int j ; for (j=1 ; j<=NX ; j++) value_assign(AK(id,j),VALUE_ZERO) ; /*j0*/ value_assign(AK(id,js),VALUE_MONE);value_assign(DK(id),VALUE_ZERO); } static void columntoline(Pproblem XX,int js, int id) /*****************/ { anonymctol(XX,js,id) ; G(id)=B(js) ; } //DN static void boundline(Pproblem XX,int io,int bou) /************/ //{ oppositeline(XX,io); DK(io) += bou; if (io>MX-NEGAL) permutel(XX,io,MX-NEGAL); //} static void boundline(Pproblem XX,int io,Value bou) /************/ { oppositeline(XX,io); value_addto(DK(io),bou); if (io>MX-NEGAL) permutel(XX,io,MX-NEGAL); } //DNstatic int addnonbasicbound(Pproblem XX,int jnb,int nbb,int nonfixe) static int addnonbasicbound(Pproblem XX,int jnb,Value nbb,int nonfixe) { int r; if (value_zero_p(nbb)) /*if (nonfixe)*/ { celimination (XX,jnb); return(0); } if ((r=newcut(XX))) return(r); anonymctol(XX,jnb,MX);boundline(XX,MX,nbb); if (!nonfixe) ++NEGAL; return(0); } //DNstatic int addedbound(Pproblem XX,int iv,int borne)//DN??? void or int??? static void addedbound(Pproblem XX,int iv,Value borne)//DN??? void or int??? { int iy; for (iy=NX ; iy>=1 ; iy--) if (B(iy)==iv) addnonbasicbound(XX,iy,borne,1); for (iy=IC1 ; iy<=MX-NEGAL ; iy++) if (G(iy)==iv) { newcut(XX); anonymouscopyline(XX,iy,MX); boundline(XX,MX,borne); } } //DN static int modifyconstr(Pproblem XX,int jfixe,int ifixe,int k,int ineq) static int modifyconstr(Pproblem XX,int jfixe,int ifixe,Value k,int ineq) { if (jfixe) { if (ineq) { int v; TMAX +=1; newcut(XX);v=G(MX);columntoline(XX,jfixe,MX);B(jfixe)=v;value_substract(DK(MX),k); if (ipivotage(XX,MX,jfixe)) XBUG(81); retirer(XX,MX); --NUMERO; if (VW5) wvoir(XX); } else addnonbasicbound(XX,jfixe,k,0); /*k==0: celimination(XX,jfixe);*/ } else { value_substract(DK(ifixe),k); if (ifixe0) c=r; else l= -r; if (VW4&&ineq) {fprintf(FTRACE,"increment y(%d)= ",v);fprint_Value(FTRACE,k);} if (modifyconstr(XX,c,l,k,ineq)) XBUG(90); if (VW5){if(c)fprintf(FTRACE,"(column)");fprintf(FTRACE," - variable %d",v);} return(0); } /* ================================================================== */ static int hb(struct rproblem *RR, int iv) { int j; for (j=1; j<=RR->nb; j++) if (RR->iname[RR->b[j]]==iv) return(j); return(0); } static void copyrline(struct rproblem *RR,int is, int id) /************/ { if (is != id) /* line is (source) is copied on line id) */ { int j ; RR->d[id]= RR->d[is] ; for (j=1 ; j<=RR->n ; j++) RR->a[id][j] = RR->a[is][j] ; } } static void rinverse(struct rproblem *RR,int ir) /************/ { int j ; RR->d[ir]= - RR->d[ir] ; for (j=1 ; j<=RR->n ; j++) RR->a[ir][j] = - RR->a[ir][j] ; } static int recherche(struct rproblem *RR, int varia) { int i; for (i=1 ; i<=RR->m ; i++) if (RR->g[i]==varia) return(-i); /* line */ for (i=1 ; i<=RR->n ; i++) if (RR->b[i]==varia) return(i); /* column */ return(0); } static void razcoutr(struct rproblem *RR, int mini) { int j; RR->d[0]=0; if (mini==1) RR->e[0]=1 ; else RR->e[0]= -1 ; for (j=1 ; j<=RR->nvar ; j++) RR->a[0][j]=0; } static int preparephase2(Pproblem XX,struct rproblem *RR,int rvar, int mini,int kost) { int ir; if (kost== 0) { if (rvar<=RR->n) ir= rvar; else ir=RR->n-rvar; } else ir=recherche(RR,rvar); if (ir<0) copyrline(RR,-ir,0); else { razcoutr(RR,mini); RR->a[0][ir] = -1 ; }; if (mini==1) RR->e[0]=1 ; else { RR->e[0]= -1 ; if (kost== -1) rinverse(RR,0); } if (kost== -1) { if (ir<0) { RR->inf[0]=RR->inf[-ir]; *(RR->pconor)=*(RR->pconor+RR->g[-ir]); } else { RR->inf[0]=0; *(RR->pconor)=*(RR->pconor+RR->b[ir]); }} RR->cost= kost; return(ir); } static int phase1(Pproblem XX,struct rproblem *RR) { int r; integertoreal(XX,RR,1); razcoutr(RR,1); if (MSR>=5) RR->meth=1; r=realsimplex(RR); RR->meth=0; return(r);/*?*/ } static int phase2(Pproblem XX,struct rproblem *RR,int rvar,int mini) { int rsr; if (MSR<=2) { integertoreal(XX,RR,1); razcoutr(RR,1); if (MSR==1) RR->meth=1; preparephase2(XX,RR,rvar,mini,0); } else { if (MSR<=3) rsr=phase1(XX,RR); /* repetee pour verifications */ if (preparephase2(XX,RR,rvar,mini,-1)==0) XBUG(51);/*introuvable*/ } if ((rsr=realsimplex(RR))) if (VIS2) messrsimplex(XX,RR,rsr,2); return(rsr); } static int dealtvariable(Pproblem XX,int v) /*********************/ { return((v<=NV && v>0)); } static int inevariable(Pproblem XX,int v) /*********************/ { return(v>NV+1 && vilbound[iv])) */ /* {symbold(XX,iv);fprintf(FTRACE,"(");fprint_Value(FTRACE,XX->ilbound[iv]);fprintf(FTRACE,") "); */ /* } */ /* fprintf(FTRACE,"\\\\\n"); */ /* } */ static int majlowbounds(Pproblem XX,struct rproblem *RR,int ope) { int i,j,iv; if (MSR<6) return(0); if (ope==0) { if(REDON)for(iv=1;iv<=NUMAX;iv++){ value_assign(XX->ilbound[iv],VALUE_MONE);value_assign(XX->ibound[iv],100000);XX->utilb[iv]=-1;} else for (iv=1;iv<=NUMAX;iv++) if(dealtvariable(XX,iv)) { value_assign(XX->ilbound[iv],VALUE_MONE); XX->utilb[iv]=-1; } /*6mai1999*/ } else if (MSR<7) { for (i=IC1 ; i<=MX ; i++) if (dealtvariable(XX,iv=G(i))) if (recherche(RR,RR->nvar+i+1-IC1)>0) value_assign(XX->ilbound[iv],VALUE_ZERO); for (j=1 ; j<=NX ; j++) if (dealtvariable(XX,iv=B(j))) if (recherche(RR,j)>0) value_assign(XX->ilbound[iv],VALUE_ZERO); } else { for (j=1 ; j<=RR->nvar ; j++) if (dealtvariable(XX,iv=RR->iname[RR->b[j]])||REDON) { if(inevariable(XX,iv)) { if (value_pos_p(XX->ilbound[iv])) XBUG(85); else if(value_neg_p(XX->ilbound[iv])&&VW5) fprintf(FTRACE,"%d:HB ",iv); } value_assign(XX->ilbound[iv],VALUE_ZERO); } /*etatfonctions(XX);*/ if (MSR>=9&&XX->state>0) /*6mai1999*/ for (i=1; i<=RR->mcontr; i++) if (dealtvariable(XX,iv=RR->iname[RR->g[i]])) if (XX->utilb[iv]) if (RR->x[RR->g[i]]>=XX->ibound[iv]) //comparison of float and Value. DN { XX->utilb[iv]=0; if(VW3) {fprintf(FTRACE,"y%d ancien ",iv);fprint_Value(FTRACE,XX->ibound[iv]); fprintf(FTRACE," actuel %13.6f\\\\\n",RR->x[RR->g[i]]);} } } return(0); } /* static int cutbounds(Pproblem XX,struct rproblem *RR,int iv, int rvar,float epss) */ /* { */ /* float rrv; int r; */ /* if (VW6) {fprintf(FTRACE,"bornes cut "); symbold(XX,iv);} */ /* if ((r=phase2(XX,RR,rvar,1))) return(r); */ /* rrv=RR->x[0]+epss; fprintf(FTRACE,"var [%d] max=%13.6f ",iv,rrv); */ /* if (rrv<1) fprintf(FTRACE," {*}CTILT {*}\\\\\n"); */ /* if ((r=phase2(XX,RR,rvar,-1))) return(r); */ /* rrv=RR->x[0]-epss; fprintf(FTRACE," min=%13.6f\\\\\n",rrv); return(0); */ /* } */ static int upperbound(Pproblem XX,struct rproblem *RR,int iv,int rvar, float epss) { float rrv; int r; Value dn_tmp; if (VW6) symbold(XX,iv); if(MSR>=7) if ((r=majlowbounds(XX,RR,1))) return(r); if (MSR>=9&&XX->state>0) { if(!XX->utilb[iv]) {if (VW3) fprintf(FTRACE,"upper inutile y%d\\\\\n",iv); return(0);} else { XX->utilb[iv]=0;if (VW3) fprintf(FTRACE,"upper effectif y%d\\\\\n",iv); }} if ((r=phase2(XX,RR,rvar,1))) return(r); if (VW6) messrsimplex(XX,RR,r,2); rrv=RR->x[0]+epss; XX->rbound[iv]=rrv;//float to Value.DN //r=XX->ibound[iv]; XX->ibound[iv]=dessous(rrv); // XX->decrement += (r-XX->ibound[iv]); //r is used just like a temporary variable, replaced by dn_tmp value_assign(dn_tmp,XX->ibound[iv]); value_assign(XX->ibound[iv],dessous(rrv)); value_addto(XX->decrement,value_minus(dn_tmp,XX->ibound[iv])); if(VW6) { if (value_lt(XX->ibound[iv],VALUE_ONE)) fprintf(FTRACE," {*} TILT {*}\\\\\n"); else {fprintf(FTRACE," bound = ");fprint_Value(FTRACE,XX->ibound[iv]);fprintf(FTRACE,"\\\\\n");}} return(0); } static int lowbound(Pproblem XX,struct rproblem *RR,int iv,int rvar, int *first,float epss) { Value range,r=VALUE_ZERO; float rrv; if(VW2&&inevariable(XX,iv)) fprintf(FTRACE,"examen fonc y%d\\\\\n",iv); if (MSR<6 || XX->ilbound[iv]) { if (MSR>=7) if ((r=majlowbounds(XX,RR,1))) return(r);; /* remplacer par 7 */ if ((r=phase2(XX,RR,rvar,-1))) return(r); rrv=RR->x[0]-epss; XX->rlbound[iv]=rrv; value_assign(XX->ilbound[iv],dessus(rrv)); value_addto(XX->decrement,XX->ilbound[iv]); if(VW2&&inevariable(XX,iv)) fprintf(FTRACE,"calcul fonc y%d\\\\\n",iv); if(VW6||(VW2&&inevariable(XX,iv))) if (MSR>=6 && (value_zero_p(XX->ilbound[iv]))) {symbold(XX,iv);fprintf(FTRACE," preuve basse par calcul\\\\\n"); }; } rrv=XX->rbound[iv]; value_assign(range,value_minus(XX->ibound[iv],XX->ilbound[iv])); if(VW6) if (value_notzero_p(XX->ilbound[iv])){ symbold(XX,iv); messrsimplex(XX,RR,r,2); fprintf(FTRACE," lbound = ");fprint_Value(FTRACE,XX->ilbound[iv]);fprintf(FTRACE," bound = "); fprint_Value(FTRACE,XX->ibound[iv]);fprintf(FTRACE," RANGE ");fprint_Value(FTRACE,range); if (value_zero_p(range)) {fprintf(FTRACE," {*} TILT {*} ");fprint_Value(FTRACE,XX->ibound[iv]);fprintf(FTRACE,"\\\\\n");} else if (value_neg_p(range)) fprintf(FTRACE," $->$ {*}{*} EMPTY {*}{*}\\\\\n"); else fprintf(FTRACE,"\\\\\n"); } if (value_neg_p(range)) return(VRVID); if (value_zero_p(range)) ++XX->tilt; if (dealtvariable(XX,iv)) if (*first||rrvrrbest) { XX->rrbest=rrv; XX->vbest=iv; *first=0;} return(0); } static Value boundssum(Pproblem XX) { int i,j,v; Value s = VALUE_ZERO; for (i=IC1 ; i<=MX ; i++) if (dealtvariable(XX,v=G(i))) { value_addto(s,XX->ibound[v]); value_substract(s,XX->ilbound[v]); } for (j=1 ; j<=NX ; j++) if (dealtvariable(XX,v=B(j))) { value_addto(s,XX->ibound[v]); value_substract(s,XX->ilbound[v]); } return(s); } static int computebounds(Pproblem XX,struct rproblem *RR,int up) { int i0,i,j,r,first; float epss; Value s1=VALUE_ZERO,s2=VALUE_ZERO; epss=0.0001; /*epss=0.000001;*/ if ((r=phase1(XX,RR))) { if(r==VRVID&&VW4) fprintf(FTRACE,"immediate "); return(r); } if (VW2) value_assign(s1,boundssum(XX)); value_assign(XX->decrement,VALUE_ZERO);XX->tilt=0;first=1;i0=RR->nvar+1-IC1; majlowbounds(XX,RR,0); if (up) { for (i=IC1 ; i<=MX ; i++) if (dealtvariable(XX,G(i))) if ((r=upperbound(XX,RR,G(i),i0+i,epss))) return(r); for (j=1 ; j<=NX ; j++) if (dealtvariable(XX,B(j))) if ((r=upperbound(XX,RR,B(j),j,epss))) return(r); } for (i=IC1 ; i<=MX ; i++) if (dealtvariable(XX,G(i))) if ((r=lowbound(XX,RR,G(i),i0+i,&first,epss))) return(r); for (j=1 ; j<=NX ; j++) if (dealtvariable(XX,B(j))) if ((r=lowbound(XX,RR,B(j),j,&first,epss))) return(r); if (REDON) { for (i=IC1;i<=MX;i++) if(inevariable(XX,G(i))) if ((r=lowbound(XX,RR,G(i),i0+i,&first,epss))) return(r); } if (!XX->tilt) { value_assign(XX->bestlow,XX->ilbound[XX->vbest]); value_assign(XX->bestup,XX->ibound[XX->vbest]);} if (VW2) value_assign(s2,boundssum(XX)); if (VW3) {fprintf(FTRACE,"etat=%d s1=",XX->state);fprint_Value(FTRACE,s1);fprintf(FTRACE," s2="); fprint_Value(FTRACE,s2);fprintf(FTRACE," diff=");fprint_Value(FTRACE,value_minus(s1,s2)); fprintf(FTRACE," decrement=");fprint_Value(FTRACE,XX->decrement);fprintf(FTRACE,"\\\\\n");} if (PMET2>=66) if (value_zero_p(XX->decrement)) XBUG(93); XX->state=1; return(0); /*xs*/ } static void voirbornes(Pproblem XX) { int iv; for (iv=1; iv<=NUMAX ; iv++) if (dealtvariable(XX,iv)) { if (presence(XX,iv)) { symbold(XX,iv); fprintf(FTRACE,"("); if (value_notzero_p(XX->ilbound[iv])) {fprint_Value(FTRACE,XX->ilbound[iv]);fprintf(FTRACE," - ");} fprint_Value(FTRACE,XX->ibound[iv]);fprintf(FTRACE,") "); } else fprintf(FTRACE," .... ");} if (REDON){ fprintf(FTRACE,"redondances: "); for (iv=1; iv<=NUMAX ; iv++) if (inevariable(XX,iv)) if (presence(XX,iv)) if (value_notzero_p(XX->ilbound[iv])) { symbold(XX,iv); fprintf(FTRACE,"(");fprint_Value(FTRACE,XX->ilbound[iv]);fprintf(FTRACE,") - "); } } fprintf(FTRACE,"\\\\\n"); /*retour*/ } static void listehb(Pproblem XX,struct rproblem *RR) { int iv,c; for (iv=1; iv<=NUMAX ; iv++) if (dealtvariable(XX,iv)) if ((c=hb(RR,iv))) if (value_le(XX->ibound[iv],VALUE_ONE)) {fprintf(FTRACE," (%d):",c); symbold(XX,iv);} fprintf(FTRACE," hors-base possibles\\\\\n"); } static int choose(Pproblem XX,struct rproblem *RR) { int i,j,iv,prem; Value bes=VALUE_ZERO; if (XX->state<0) XBUG(61);if (XX->state==0) XBUG(62); prem=1;XX->tilt=0; for (i=IC1 ; i<=MX ; i++) if (dealtvariable(XX,iv=G(i))) if (prem||(value_lt(XX->ibound[iv],bes))) { value_assign(bes,XX->ibound[iv]); XX->vbest=iv; prem=0;}//bes initialized by prem = true :-( for (j=1; j<=NX; j++) if (dealtvariable(XX,iv=B(j))) if (prem||(value_lt(XX->ibound[iv],bes))||((MSR>=8)&&(RR->state>=0)&&(value_eq(XX->ibound[iv],bes))&&(hb(RR,iv)))) {value_assign(bes,XX->ibound[iv]);XX->vbest=iv;prem=0;} value_assign(XX->bestlow,XX->ilbound[XX->vbest]); value_assign(XX->bestup,XX->ibound[XX->vbest]); if (prem) XBUG(63); if (VW4&&(value_gt(XX->bestup,int_to_value(PCOMP)))) /*if (VIS2&&XX->bestup>=2)*/ { //DN. int_to_value symbold(XX,XX->vbest);fprintf(FTRACE,"DILATATION = ");fprint_Value(FTRACE,XX->bestup); voirbornes(XX);} return(0); } static int ajout(Pproblem XX) { int iv,i; Value bi; if (REDON) for (i=MX;i>=IC1;i--) if (inevariable(XX,iv=G(i))) if (XX->ilbound[iv]) { if(VW4) fprintf(FTRACE,"retrait inequa %d\\\\\n",iv); retirerineq(XX,i); } if (PMET2<7) return(0); razcut(XX); for (iv=1 ; iv<=NUMAX ; iv++) { if (dealtvariable(XX,iv)) { if (value_pos_p(value_assign(bi,int_to_value(XX->ibound[iv])))) addedbound(XX,iv,bi); else { if (value_neg_p(bi)) XBUG(72); else {if(VW6) symbold(XX,iv);} }} } if(VW7) wvoir(XX); return(0); } static int majuu(Pproblem XX,Pproblem UU,struct rproblem *RR) { int iv,ns; Value bl; ns=0; for (iv=1; iv<=NUMAX; iv++) if (dealtvariable(XX,iv)) { if (value_pos_p(value_assign(bl,XX->ilbound[iv]))) { if (VW4||(VW3&&MSR>=8)) {fprintf(FTRACE,"majuu, y%d=",iv);fprint_Value(FTRACE,bl);fprintf(FTRACE,"\\\\\n");} if (findandmod(UU,iv,bl,1)) XBUG(95); vidr(RR); } value_assign(UU->ibound[iv],value_minus(XX->ibound[iv],XX->ilbound[iv])); value_assign(UU->ilbound[iv],VALUE_ZERO); if (value_ge(UU->ibound[iv],2)) ++ns; } return(0); } static int majb(Pproblem XX,Pproblem UU,struct rproblem *RR) { if (PMET2>=3) {/*2MMM3*/ if (majuu(XX,UU,RR)) XBUG(97); UU->state=1; /*xs*/ if(VW4) { fprintf(FTRACE,"apres majuu: "); voirbornes(UU);} if (satisfait(UU)) { if (VW4) wvoir(UU); return(99); } /*(VRFIN)*/ } else { copybound(XX,UU); UU->state=1;} /*xs*/ return(0); } static int failuresp(Pproblem XX,struct rproblem *RR,int *compb) { int r,up; *compb=1; XX->bloquer=0; up=1; if (VW7) {fprintf(FTRACE," failure:\n"); wvoir(XX); } if (PMET2>=66) XX->bloquer=1; else if (PCOMP&&XX->state>0) { if ((r=choose(XX,RR))) return(r); if (XX->tilt) XBUG(67); if ((PCOMP==9)||(value_le(XX->bestup,VALUE_ONE))) *compb=0; /*XX->bestup<=PCOMP*/ else if (PCOMP==2) up=0; } if (!(*compb)) { if (MSR>=8&&RR->state==0&&RR->vc){ if (VW3) fprintf(FTRACE,"SIMULER y%d col%d\\\\\n",RR->namev,RR->vc); if (!RR->vc) XBUG(84); if (RR->vc!=hb(RR,RR->namev)) XBUG(87); elimination(RR,RR->vc); } else if ((r=phase1(XX,RR))) { if (r==VRVID&&VW4) fprintf(FTRACE,"qimmediate "); return(r);} if ((r=choose(XX,RR))) return(r); XX->state++; /*xs*/ if (VW3&&MSR>=8) { if (hb(RR,XX->vbest)) fprintf(FTRACE,"HB:%d\\\\\n",XX->vbest); else fprintf(FTRACE,"echec HB:%d\\\\\n",XX->vbest);} } else { if (VW3&&MSR>=8) fprintf(FTRACE,"COMPUTE TOTAL %d\\\\\n",XX->state); if ((r=computebounds(XX,RR,up))) return(r); if (VW4) {fprintf(FTRACE,"apres compute:"); voirbornes(XX);} } return(0); } static int reducedpb(Pproblem XX,Pproblem VV,struct rproblem *RR) { ++VV->repere; if (VW2) { int tete,i; if (NV>10) tete=NV-NX; else tete=10-NX; for (i=1; i<=tete; i++) fprintf(FTRACE,"---"); fprintf(FTRACE,"(%d Variables) REDUCED PROBLEM %d equa ",NX,NEGAL); for (i=MX+1-NEGAL; i<=MX; i++) symbold(XX,G(i)); if (VW4) fprintf(FTRACE," repere %d\\\\\n",VV->repere); else fprintf(FTRACE,"\\\\\n"); if (VW5) wvoir(XX); } ICOUT=0; VRESULT=is2(XX,VV,RR,1) ; if (VIS2) if ((VRESULT!=VRFIN) && (VRESULT!=VRVID)) fprintf(FTRACE,"DEAD reduced pb NX=%d NEGAL=%d vr=%d\\\\\n",NX,NEGAL,VRESULT); return(VRESULT) ; } static int splitpb(Pproblem XX,Pproblem VV,struct rproblem *RR, int jfixe,int ifixe, Value k,int ineq) { int vhb; ++VV->repere; if (VW2) { int tete,i,va; if (NV>10) tete = NV - NX; else tete = 10 - NX; for (i=1; i<=tete; i++) fprintf(FTRACE,"---"); fprintf(FTRACE,"{.} (%d Variables) SPLIT PROBLEM ",NX); if (VW4) fprintf(FTRACE," repere %d ",VV->repere); if (jfixe) va = B(jfixe) ; else va = G(ifixe); fprintf(FTRACE,"$");symbol(XX,va); if (ineq) {fprintf(FTRACE,", cost <= ");fprint_Value(FTRACE,value_uminus(k));fprintf(FTRACE,"$\\\\\n");} else fprintf(FTRACE," = ");fprint_Value(FTRACE,k);fprintf(FTRACE,"$\\\\\n"); } if (ineq) { copyline(XX,1,++MX); if (modifyconstr(XX,jfixe,MX,k,ineq)) XBUG(91); } else if (modifyconstr(XX,jfixe,ifixe,k,ineq)) XBUG(92); if (VW4&&(ineq||VW6)) wvoir(XX); if (MSR>=8&&RR->state==0) { if ((vhb=hb(RR,XX->vbest))) { if (RR->vc!=hb(RR,RR->namev)) XBUG(88);if (XX->vbest!=RR->namev) XBUG(89); } else if (VW3) fprintf(FTRACE,"NON, car basic:%d\\\\\n",XX->vbest); } ICOUT=0; VRESULT=is2(XX,VV,RR,1) ; if (VIS2) if (VRESULT!=VRFIN && VRESULT!=VRVID) { fprintf(FTRACE,"DEAD split pb NX=%d vr=%d k=",NX,VRESULT);fprint_Value(FTRACE,k);fprintf(FTRACE,"\\\\\n");} return(VRESULT) ; } static int recup(Pproblem XX,Pproblem VV) /*2MMM4*/ { if (PMET2>=4) if (VV->marque==0) { if (copystructure(XX,VV)) XBUG(43); VV->marque=1;} return(0); } static int reduction(Pproblem XX,struct rproblem *RR) { int iv; Value k; vidr(RR); for (iv=1; iv<= XX->numax; iv++) if (dealtvariable(XX,iv)) if (value_eq((value_assign(k,XX->ibound[iv])),XX->ilbound[iv])) if (findandmod(XX,iv,k,0)) XBUG(98); return(ajout(XX)); } static int eclater(Pproblem XX,Pproblem UU,Pproblem VV, struct rproblem *RR) { int r,ivb; Value k; struct problem AYY; vid(&AYY); /*ssss*/ if ((r=recup(UU,VV))) return(r); if ((r=choose(UU,RR))) return(r); ivb=UU->vbest; if (value_notzero_p(UU->ilbound[ivb])) XBUG(54); for (value_assign(k,UU->ilbound[ivb]); value_le(k,UU->ibound[ivb]); value_increment(k)) { if (copystructure(UU,&AYY)) XBUG(44); AYY.vbest=ivb; if ((r=ajout(&AYY))) return(r); if (value_notzero_p(k)) vidr(RR); RR->namev=ivb; RR->vc = hb(RR,ivb); if (!presence(&AYY,ivb)) XBUG(53); r=splitpb(&AYY,VV,RR,AYY.cu,AYY.lu,k,0); if (r==VRFIN) value_assign(XX->dk[1],AYY.dk[1]); if (r!=VRVID) return(r); } if (VW4) fprintf(FTRACE,"INDIRECT PROOF NX=%d\\\\\n",NX); return(VRVID); } static int failure(Pproblem XX,Pproblem UU,Pproblem VV, struct rproblem *RR) { int r,compb; if ((r=failuresp(XX,RR,&compb))) { if (r!=VRVID) XBUG(99); if (VW4) fprintf(FTRACE,"real proof\\\\\n"); return(r);} if (compb) if ((r=majb(XX,UU,RR))) { if (r!=99) XBUG(96); else { value_assign(XX->dk[1],UU->dk[1]); if ((r=recup(UU,VV))) XBUG(94); if (VW4) fprintf(FTRACE,"integ sol\\\\\n"); return(VRFIN); }} if (VW3&&MSR>=8) {fprintf(FTRACE,"U: "); listehb(UU,RR);} if (XX->tilt||XX->bloquer) { if ((r=reduction(UU,RR))) return(r); if ((r=reducedpb(UU,VV,RR))==VRFIN) { if (VW4&&XX->bloquer) fprintf(FTRACE,"solution bloquee\\\\\n"); value_assign(XX->dk[1],UU->dk[1]); } return(r); } else return(eclater(XX,UU,VV,RR)); } static int fastplus(Pproblem XX,Pproblem VV,struct rproblem *RR) { int r; struct problem AUU; vid(&AUU); /*ssss*/ if (PMETH>=8) return(failure(XX,XX,VV,RR)); /*PPP8 mise au point*/ if (copystructure(XX,&AUU)) XBUG(45); /*revoir utilite*/ if (PMETH>=7) return(failure(XX,&AUU,VV,RR)); /*PPP7*/ NITER=0; if (CRITERMAX) TMAX=(NX * CRITERMAX)+NITER; if (PMETH>=6) /*PPP6*/ { if ((r=dualentier(XX))==VRINS) return(failure(XX,&AUU,VV,RR)); if (VW2) if (r==VRFIN){ fprintf(FTRACE,"Feasibility proof (dual), cost = "); fprint_Value(FTRACE,value_uminus(DK(1)));fprintf(FTRACE,"\n"); vzlast(XX); if (VW4) wvoir(XX); } } else { struct problem AZZ; vid(&AZZ); /*ssss*/ if ((r=fast(XX,&AZZ))==VRINS) { if (PMET3==0) r=failure(XX,&AUU,VV,RR); else { r=failure(&AZZ,&AUU,VV,RR); value_assign(XX->dk[1],AZZ.dk[1]); } ICOUT=0; /******* revoir!!!!! necessaire 2 cas ********/ } } return(r); } static int pseudocopie(Pproblem XX,Pproblem YY,int condition) { Pproblem pt; if (XX->state<0) XBUG(42); if (condition) {pt=XX;XX=YY;YY=pt;} else if (copystructure(XX,YY)) XBUG(47); return(0); } static int iprimalplus(Pproblem XX,struct rproblem *RR,int minimum,int i2,int stopcout) { int r,step,dispo;float rsr; Value bornesr,bornepai,bornea,solution; struct problem AYY,AZZ,AVV,AWW,*pty,*ptv,*ptw; vid(&AYY); vid(&AZZ); vid(&AVV); vid(&AWW); dispo=0; /*ssss*/ /*xs*/ pty=&AYY;ptv=&AVV;ptw=&AWW; XX->state=0;if (copystructure(XX,&AZZ)) XBUG(46); if ((r=iprimal(XX,minimum,i2,stopcout))!=VRINS) return(r); value_assign(bornepai,DK(ICOUT)); rsr=localrsimplex(XX,1); if (XX->minimum) value_assign(bornesr,dessous(rsr)); else value_assign(bornesr,dessus(rsr)); if (VW2) {fprintf(FTRACE,"PRIMAL FAILS (%d ITER), possible integer values %13.6f $<$ ",NITER,-rsr); fprint_Value(FTRACE,value_uminus(bornesr));fprintf(FTRACE," $<=$ cost $<=$ "); fprint_Value(FTRACE,value_uminus(DK(ICOUT)));fprintf(FTRACE,"\\\\\n");} value_assign(solution,bornepai); value_assign(bornea,bornesr); VRESULT=VRFIN; step=0; while (value_lt(solution,bornea)) { Value possible,expossible=VALUE_ZERO,typspl; ++step; /*if (XX->minimum)*/ if (value_eq(value_plus(solution,VALUE_ONE),bornea)) value_assign(possible,value_plus(solution,VALUE_ONE)); else if (PDICO==1 || (step==1&&PDICO==2)) value_assign(possible,value_div(value_plus(bornea,solution),2)); else if (PDICO==3) value_assign(possible,value_plus(solution,VALUE_ONE)); /* descendant */ else if (PDICO==2 && VRESULT==VRFIN) value_assign(possible,value_plus(solution,VALUE_ONE)); else value_assign(possible,bornea); /* ascendant. Par curiosite, si 5, ->inequation */ if ((value_eq(possible,bornea)) && (PDICO!=5)) typspl=0; else typspl=1; if (VW2) {fprintf(FTRACE,"dichoto%d NX=%d poss=",step,NX);fprint_Value(FTRACE,possible);fprintf(FTRACE,"\\\\\n");} ptw->marque=0; ptw->repere=11; if (dispo) { /*2MMM5*/ if (pseudocopie(ptv,pty,(PMET2>=5 && (value_eq(possible,value_plus(solution,VALUE_ONE)))))) XBUG(47); findandmod(pty,G(1),(value_minus(possible,expossible)),1); /*provisoire*/ if (VW4) wvoir(pty); vidr(RR); VRESULT=reducedpb(pty,ptw,RR); } else { if (copystructure(&AZZ,pty)) XBUG(48); vidr(RR); VRESULT=splitpb(pty,ptw,RR,0,ICOUT,possible,typspl); } if (VRESULT!=VRFIN&&VRESULT!=VRVID) return(VRESULT); if (PMET2>=4) { if(PMET2<7)if (!ptw->marque) XBUG(55);if (VIS2&&PMET2<7) wvoir(ptw);} if (VRESULT==VRVID) value_assign(bornea,value_minus(possible,VALUE_ONE)); else { if (pseudocopie(ptw,ptv,(PMET2>=5 && dispo))) XBUG(49); /*2MMM5*/ if (PMET2>=4) {dispo=1; value_assign(expossible,possible); } /*2MMM4*/ if (typspl==0) value_assign(solution,possible); else if (value_eq(value_assign(solution,pty->dk[1]),possible)) { if (VIS2) { fprintf(FTRACE," p = ");fprint_Value(FTRACE,value_uminus(possible)); fprintf(FTRACE," s = ");fprint_Value(FTRACE,value_uminus(solution));} XBUG(56);} else if (VW1) { fprintf(FTRACE," verified solution ");fprint_Value(FTRACE,value_uminus(possible)); fprintf(FTRACE," improved solution ");fprint_Value(FTRACE,value_uminus(solution)); fprintf(FTRACE," \\\\\n"); } } /*if(PMETH==8) XBUG(39);*//*mesure cout step 1*/ } if (VW1) {fprintf(FTRACE,"Final solution:");fprint_Value(FTRACE,value_uminus(solution));fprintf(FTRACE,"\\\\\n");} return(VRESULT); } /* ======================================================================= */ /* DYNAM.c */ /* ======================================================================= */ int dynam(Pproblem XX,int nd) /****************************/ { /*struct problem AVV; struct rproblem ARR; vid(&AVV); vidr(&ARR);*/ /*fprintf(FTRACE,"fourier=%d varc=%d choixpiv=%d forcer=%d\\\\\n", PFOURIER,VARC,CHOIXPIV,FORCER);*/ int ca,v,uni; ++XX->dyit; uni=0; if (DYN==1) { printf("dynamic visualization\n"); XX->itdirect=0; DYN=2; return(0); } if (XX->niter < XX->itdirect) return(0); printf("passage dynamic iter no=%d ----------------\n",XX->dyit); if (nd==1) printf("dual dual dual dual dual dual iteration\n"); if (DYN) printf("dynam, iteration %d ",NITER); printf("query: "); ca= getchar(); fclose(XX->ftrace) ; while (ca!=10) { if (ca==117) { uni=1;printf("unimodular change "); } else if (ca==118) { printf("visu= "); scanf("%2d",&v); VISU=v;} else if (ca==104) /*h*/ { printf("help\n"); printf(" visu=%d\n",VISU); printf(" direct "); if(DYN)printf("NO\n");else printf("YES\n"); printf(" itdirect=%d\n",XX->itdirect); } else if (ca==100) DYN=0; /*d*/ else if (ca==115) /*s*/ { printf("System state will be vizualized\n"); XX->ftrace=fopen("/tmp/jtrace.tex","a") ; wvoir(XX) ; fclose(XX->ftrace) ; } else if (ca==105) /*i*/ { printf("go to iteration: "); scanf("%2d",&v); XX->itdirect=v; DYN=2;} else if (ca!=10) printf("commande %d?\n",ca); /*else if (ca!=10) printf("commande ?\n");*/ printf("another query: "); ca= getchar(); ca= getchar(); } XX->ftrace=fopen("/tmp/jtrace.tex","a") ; return(uni); /************************* while (c!=10) { c= getchar(); printf("vidage caractere =%d\n",c); }************/ /*printf("selon c c=%c\n",c); printf("selon d c=%d\n",c);*/ /*v=isalpha(c); printf("produit par isalpha v=%d\n",v);*/ /*v=isdigit(c); printf("produit par isdigit v=%d\n",v);*/ /*printf("conversion \n"); v=atoi(&c); printf("produit par atoi v=%d\n",v);*/ /*printf("action direct:0 suite:1 visu:2 dynam:3 ? ");scanf("%2d",&v);*/ } /* ======================================================================= */ /* RESULT PROCEDURE */ /* ======================================================================= */ static int vzlast(Pproblem XX) /********************************/ { if (VISU>=ZVL2) wvoir(XX); if (VISU>=ZVL2) fprintf(FTRACE,"******** "); if (VRESULT==VRVID) fprintf(FTRACE,"empty polyhedron"); else if (VRESULT==VRBUG) fprintf(FTRACE,"bug programme %3d",NUB); else if (VRESULT==VRINF) fprintf(FTRACE,"solution infinie"); else if (VRESULT==VRFIN) fprintf(FTRACE,"integer solution"); else if (VRESULT==VRDEB) { fprintf(FTRACE,"debordement tableaux "); if (MX >MAXMX) fprintf(FTRACE," line MX=%4d\n",MX); else fprintf(FTRACE," column NX=%4d\n",NX); } else if (VRESULT==VRINS) fprintf(FTRACE,"trop de pivotages"); else if (VRESULT==VROVF) fprintf(FTRACE,"overflow"); else if (VRESULT==VRCAL) fprintf(FTRACE,"appel errone **"); else fprintf(FTRACE,"bug message * result=%5d",VRESULT); fprintf(FTRACE," after %d iterations\\\\\n",NITER) ; return 0; } /* ======================================================================= */ /* MAIN PROCEDURE */ /* ======================================================================= */ static int is2(Pproblem XX,Pproblem VV,struct rproblem *RR) { int nlibre; /* nlibre: information donnee par step 4 au step 5 */ int prevnx,stfou; /*int forcer; forcer=(PFOURIER>3);*/ /* ICOUT: ligne ou se trouve la fonction cout a optimiser; information indiquee par le step 1 (initialisations); en absence de cout, =0 */ /* 111111111111111111111111 VERIFICATIONS INITIALISATIONS 11111111111111 */ /* 22222222222222222222 ELIMINATION EGALITES 22222222222222222222222222 */ nlibre =0; prevnx = 0; stfou = 0; if (VISU) if (NEGAL) vzstep(XX,1); while (NEGAL) { int r,ik=0,jk; Value pivot; { int i; Value dm = VALUE_MONE; for (i=MX ; i>=1 ; i--) if (janus_egalite(XX,i)) { int j; Value sk,de; value_assign(pivot,VALUE_ZERO); for (j=NX; j>=1; j--) { if (value_assign(sk,value_abs(AK(i,j)))) if (value_one_p(value_assign(pivot,sk))) { ik=i; jk=j; break;}} if (value_zero_p(pivot)) /* equation surabondante ou incompatible */ { if (value_notzero_p(DK(i))) return(VRVID); if (VISU>=ZVEV) fprintf(FTRACE, "Equation $e_{%d}$ is empty.\\\\\n",G(i)); retirer(XX,i); NEGAL-- ; break ; /*t603 */ } /* next line : t602 */ if (value_one_p(pivot)) break; if ((value_lt(value_assign(de,value_abs(DK(i))),dm))||(value_neg_p(dm))) { value_assign(dm,de); ik=i; /*choice of the equality with the smallest rhs*/ } } } if (value_zero_p(pivot)) continue ; if (value_notone_p(pivot)) { /* unimodular change of variables is necessary */ int jojo; if (gcdtest(XX,ik,&jk,&pivot,ZVG2)) return(VRVID); value_assign(pivot,VALUE_ZERO); for (jojo=NX; jojo>=1; jojo--) if (value_one_p(value_abs(AK(ik,jojo)))) { jk=jojo; value_assign(pivot,VALUE_ONE); break;} if (value_notone_p(pivot)) if ((r=build0(XX,ik,1,&jk))) return(r); } else if (VISU>=ZVEV+1) fprintf(FTRACE,"Equation $e_{%d}$ includes a unitary coefficient.\\\\\n",G(ik)); if (VISU>=ZVEV+1) igvoir3(XX,1,MX,1,ik,jk,0,0,0,0,0); if ((r=ipivotage(XX,ik,jk))) return(r); celimination(XX,jk); if (fluidline(XX,ik)) retirer(XX,ik); else permutel(XX,ik,MX+1-NEGAL); NEGAL-- ; }//of while if (NEGAL) XBUG(15); /* equalities still remain */ /* 2323232323232323232323 INEQUALITIES DIVISION 23232323232323232323 */ { int i; if (VISU) vznewstep(XX,2); for (i=IC1 ; i<=MX ; i++) { inequaldivision(XX,i); } } /* 3333333333333333333333 FOURIER-MOTZKIN trivial 33333333333333333333333 */ if (PFOURIER) { int i,j; if (VISU) vznewstep(XX,3); for (j=NX ; j>=1 ; j--) if (freecolumn(XX,j)) if (possible(XX,j)) { int tplus,tmoins; Value k; tplus= tmoins= 0 ; for (i=MX ; i>=IC1 ; i--) if (value_assign(k,AK(i,j))) {if (value_pos_p(k)) tplus++ ; else tmoins++ ;} if (tplus==0||tmoins==0) { if (tplus+tmoins) { fourier0(XX,j); if (!FORCER) if (ICOUT==0 && satisfait(XX)) return (VRFIN); } else { /* ********** colonne vide */ if(VISU>=ZVFEC) vzemptycol(XX,j); emptycolelim(XX,j); } } } } /* 4444444444444444444 FOURIER-MOTZKIN - DUMMY ELIMINATIONS 4444444444444*/ stfou=0;// fprintf(stderr,"%d %d %d\n",nlibre,prevnx,stfou); while (stfou==0 || (NX=1 ; j--) if (freecolumn(XX,j)) { int tplus,tmoins,ip1=0,im1=0,tp2,tm2,i1; Value sk,dkmin; tplus= tmoins= tp2= tm2=0 ; i1= -100000 ; value_assign(sk,VALUE_ZERO); value_assign(dkmin,int_to_value(1000000)); for (i=MX ; i>=IC1 ; i--) if (value_notzero_p(AK(i,j))) { value_assign(sk,AK(i,j)); if (value_pos_p(sk)) { tplus++ ; {if (value_one_p(sk)) { ip1= i ; if (value_lt(DK(i),dkmin)) { i1=i; value_assign(dkmin,DK(i)); } } else tp2++;} } else { tmoins++ ; {if (value_mone_p(sk)) { im1=i;if (value_lt(DK(i),dkmin)) {i1=i; value_assign(dkmin,DK(i)); } } else tm2++ ;}} } if (value_zero_p(sk)) /* ********** colonne vide */ { if (possible(XX,j)) { if (VISU>=ZVFEC) vzemptycol(XX,j); emptycolelim(XX,j); } else if (VISU>=ZVFW) fprintf(FTRACE,"****colonne%3d vide bloquee\n",j);/*provisoire*/ } else { if (i1 < 0) { nlibre++ ; continue ; } /* ***** pas de pivot unitaire possible */ if (possible(XX,j)) { if (!FORCER) if (ICOUT==0 && satisfait(XX)) return (VRFIN) ; if (PFOURIER>0) if (tplus==0||tmoins==0) { fourier0(XX,j); continue; } if (PFOURIER>1) { if (tp2==0||tm2==0) /* criteres F-M 1 et 2 */ { if (tplus==1&&tp2==0) {fourier1(XX,1,j,ip1); continue;} else if (tmoins==1 && tm2==0 ) { fourier1(XX, -1,j,im1); continue; } else if (PFOURIER>2 && tmoins==2 && tplus==2 ) { fourier22(XX,j) ; continue ; }}} } if (VARC) /* dummy elimination */ { if (VISU>2) fprintf(FTRACE,"Dummy elimination: "); if (ipivotage(XX,i1,j)) return(VRESULT) ; retirer(XX,i1) ; if (VISU>=ZVTS) wvoir(XX); /*mauvais*/ } else nlibre++ ; }} } } /* ------------------------------- */ if (!FORCER) if (ICOUT==0 && satisfait(XX)) return (VRFIN) ; if (NX==0 && (!satisfait(XX))) return(VRVID) ; /* envisager traitement cas trivial NX==1 */ /* 555555555555555 NON-NEGATIVE VARIABLES 55555555555555 */ // fprintf(stderr,"%d %d %d\n",nlibre,prevnx,stfou);//DNDB if (!nlibre) NTP=NX; else { int i,j; if (VISU) vzstep(XX,4); if (VARC<=2) /******** strategie simplexe *******/ { int jj=0,idk=0; Value dkn = VALUE_ZERO; NTP=0; /** regroupement des variables contraintes **/ for (jj=1 ; jj<=NX ; jj++) if (!freecolumn(XX,jj)) permutec(XX,++NTP,jj); if (VISU>=ZVNPC) wvoir(XX); for (;;) { if (NTP==NX) break; if (!FORCER) if (ICOUT==0 && satisfait(XX)) return (VRFIN) ; // value_assign(dkn,VALUE_ZERO); value_assign(dkn,int_to_value(1000000)); /*provisoire*/ // fprintf(stderr,"dkn = int_to_value(1000000) = ");fprint_Value(stderr,dkn);fprintf(stderr,"\n");//DNDB for (i=MX ; i>=IC1 ; i--) if (value_lt(DK(i),dkn)) { value_assign(dkn,DK(i)); idk=i; } jj=NTP+1; if (VISU>=ZVN1) { fprintf(FTRACE,"rhs=");fprint_Value(FTRACE,dkn);fprintf(FTRACE,", let us select inequality $y_{%d}$.\\\\\n",G(idk)); } if (build1(XX,idk,jj) <0) { /* segment libre est vide! */ int absent ; absent=1; if (VISU>=ZVN1) fprintf(FTRACE,"No non-zero coefficient!!\\\\\n"); for (i=MX ; i>=IC1 ; i--) if ((absent||(value_lt(DK(i),dkn))) && jnsegment(XX,i,jj,NX)>=0) { value_assign(dkn,DK(i)); idk=i; absent=0; } if (absent) { if (VISU>=ZVNW) fprintf(FTRACE,"WARNING, all columns of free variables are empty\\\\\n"); if (ICOUT==0) NX=NTP; break; } if (VISU>=ZVN1) fprintf(FTRACE,"Then, we select inequality $y_{%d}$.\\\\\n",G(idk)); if (value_pos_p(dkn)) if (VISU>=ZVN2) fprintf(FTRACE,"WARNING SIMPLEX STRATEGY rhs positif\\\\\n"); if (build1(XX,idk,jj) <0) XBUG(33); } ++NTP; // fprintf(stderr,"%d %d %d %d %d %d %d",nlibre,prevnx,stfou, i,j,jj,idk);fprint_Value(stderr,dkn);fprintf(stderr,"\n");//DNDB if ((i=cutiteration(XX,idk,jj,0))) return(i); }//of for } else if (VARC==3) { /********* variable unique *************/ if (VISU>1) fprintf(FTRACE,"SINGLE ADDED VARIABLE\n"); if (++NX>MAXNX) return(VRDEB) ; for (i=MX ; i>=1 ; i--) { Value temporaire ; value_assign(temporaire,VALUE_ZERO); for (j=1 ; j3) wvoir(XX) ; } else if (VARC>=4) {/******* variables libres doublees *******/ for (j=NX ; j>=1 ; j--) if (freecolumn(XX,j)) { if (++NX>MAXNX) return(VRDEB); B(NX)=B(j) ; /* tres provisoire */ for (i=MX; i>=1; i--) value_assign(AK(i,NX),value_uminus(AK(i,j))); } if (VISU) fprintf(FTRACE,"map: variables libres ont ete doublees \n"); } } /* 5656565656565656 FOURIER ON CONSTRAINED VARIABLES 5656565656565656 */ if (PFOURIER) { int i,j; if (VISU) vznewstep(XX,5); for (j=NX ; j>=1 ; j--) if (possible(XX,j)) { /*if (freecolumn(XX,j))*/ int tplus,tmoins; Value k; tplus= tmoins= 0 ; for (i=MX ; i>=IC1 ; i--) if (value_assign(k,AK(i,j))) { if (value_pos_p(k)) tplus++ ; else tmoins++ ;} if (tplus==0) { if (VISU>=ZVS) vznowstep(XX); if (tplus+tmoins) { fourier0(XX,j); if(!FORCER) if (ICOUT==0 && satisfait(XX)) return(VRFIN);} else { if(VISU>=ZVFEC) vzemptycol(XX,j); emptycolelim(XX,j); } /* ********** colonne vide */ } } } /* 66666666666666 VARIOUS REDUNDANT INEQUALITIES 666666666666666 */ if (XX->remove) { int i,u,i1,i2,exmx; exmx=MX; if (VISU>=ZVS) vznewstep(XX,6); for (i=MX ; i>=IC1 ; i--) /******** empty lhs (left-hand side */ if (emptylhs(XX,i)) { if (value_neg_p(DK(i))) return(VRVID); VREDUN(++NREDUN)=G(i); if(VISU>=ZVR1) vzredundant(XX,i,-1); retirer(XX,i); } for (i=MX ; i>=IC1 ; i--) /************ intrinsic redundancies */ if (useless(XX,i)==1) {if (VISU>=ZVR1) vzredundant(XX,i,0); VREDUN(++NREDUN)=G(i);retirer(XX,i);} for (i1=MX ; i1>=IC1+1 ; i1--) /****** redundancies when comparison */ for (i2=i1-1 ; i2>=IC1 ; i2--) if ((u=redundant(XX,i1,i2))) { if (VISU>=ZVR1) { int v; if (u==i1)v=i2;else v=i1; vzredundant(XX,u,v);} VREDUN(++NREDUN)=G(u); retirer(XX,u); break; } if ((VISU>=ZVR3) && MXnredun!=exmx-MX) XBUG(69); } /* 66666666666666666 RESOLUTION SIMPLEXE ENTIERS phase 1 666666666666666666 */ { int r; if (satisfait(XX)) { r=VRFIN; if (VISU>=ZVS) fprintf(FTRACE,"{\\bf {*} PHASE: DUAL SIMPLEX} is unnecessary\\\\\n"); } else { if (PMETH==MTDAI||PMETH==MTDI2) r=dualentier(XX); else r=fastplus(XX,VV,RR); } if ((r!=VRFIN) || !ICOUT) return(r); } /* 777777777777777777 OPTIMISATION FONCTION COUT ENTIERE 777777777777777777 */ if (VW2) { fprintf(FTRACE,"Primal %d variables precr=%d comp=%d msr=%d dico=%d met3=%d met2=%d meth=%d - varc=%d choixpiv=%d crit=%d m8=%d\\\\\n",NX,PRECR,PCOMP,MSR,PDICO,PMET3,PMET2,PMETH,VARC,CHOIXPIV,CRITERMAX,REDON); wvoir(XX); } // fprintf(stderr,"%d %d %d\n",nlibre,prevnx,stfou);//DNDB return(iprimalplus(XX,RR,XX->minimum,ICOUT,0)); /* ./shj p-6 ./shj p3-7 */ } static int init_janus(Pinitpb II,Pproblem XX,int suit1) { /* ICOUT: ligne ou se trouve la fonction cout a optimiser; information indiquee par le step 1 (initialisations); en absence de cout, =0 */ /* 111111111111111111111111 VERIFICATIONS INITIALISATIONS 11111111111111 */ NITER=0; TMAX=2*MC+NV; VRESULT=0 ; XX->nredun=0; if (VISU>=ZVS) fprintf(FTRACE, "\\begin{center} {\\bf ENTER INTEGER SOLVER} \\end{center}\n"); if (MC>MAXLIGNES || NV>MAXCOLONNES) { if (VISU) fprintf(FTRACE,"entree: mc=%d maxl=%d --- nv=%d maxcol=%d\n",MC,MAXLIGNES,NV,MAXCOLONNES); return(VRDEB); } if (NP>NV) return(VRCAL) ; if (!suit1) {/* equations a la fin, inequations m sens */ int i,j,i1,i2,ifi,ni,ii2,mf,ncou,Ei; NUMAX=NUMERO=NV+MC ; NX=NV; MX=MC; mf=ni=0 ; ncou=0; ICOUT=0; LASTFREE=0;//ii2=0; for (j=1 ; j<=NV ; j++) B(j)=j ; for (i=1 ; i<=MC ; i++) { Ei = E(i); //if (ii2!=0) fprintf(stderr,"Ei = %d, E(i) = %d, I->e[i] = %d",Ei,II.e(i),II.e[i]);//DNDB if (( abs(Ei)==2)||(Ei==3)) { ++mf; //number of function "cout" DN if (( abs(Ei) ==2)&&(++ncou>1)) { if (VISU>=ZVB) fprintf(FTRACE,"nombre de couts>1\n"); return(VRCAL); } } else if (Ei!=0) { if (Ei!=1 && Ei!= -1) { if (VISU>=ZVB) fprintf(FTRACE,"! e[%3d]=%3d\n",i,Ei); return(VRCAL); } // if (ii2!=0) fprintf(stderr,"Ei danglefai+= %d",Ei);//DNDB ni++ ;//number of inequation DN } }//of for ii2=0;ifi=0;i1=mf;i2=mf+ni; IC1=mf+1; ii2=IC1-1+ni;/*derniere inequation*/ NEGAL=MX-ii2; for (i=1 ; i<=MC ; i++) { int i3; i3=(abs(Ei=E(i))==2)? mf :(Ei==3)? ++ifi :(Ei) ? ++i1 :++i2; if (Ei== -1) { value_assign(DK(i3),value_uminus(D(i))); for (j=1; j<=NV; j++) value_assign(AK(i3,j),value_uminus(A(i,j))); } else { value_assign(DK(i3),value_uminus(D(i))); for (j=1; j<=NV; j++) value_assign(AK(i3,j),value_uminus(A(i,j))); } G(i3)=NV+i ; if (abs(Ei)==2) { ICOUT=i3; XX->minimum=(Ei==2); } } if (VISU) { VDUM=0; if (VISU>=ZVA1) fprintf(FTRACE," %d variables %d functions\ %d inequalities %d equalities\\\\\n",NX,mf,ni,MC-ni-mf); if(VISU>=ZVA4) wvoir(XX); /*stepvoir(XX);*/ } //fprintf(stderr,"\n DN %d %d %d %d %d %d %d %d %d %d",i,j,i1,i2,ifi,ni,ii2,mf,ncou,Ei); //DNDB } //of if suit1 return(0); } /**********************************************************************/ int isolve(Pinitpb II, Pproblem XX,int suite) /***********/ { struct problem AVV;struct rproblem ARR;int suit1,suit2,su3,su,r,stopc; vid(&AVV);vidr(&ARR); /*fprintf(FTRACE,"fourier=%d varc=%d choixpiv=%d forcer=%d\\\\\n", PFOURIER,VARC,CHOIXPIV,FORCER);*/ if (DYN) dynam(XX,0); su3=suite/100; su=suite-100*su3; suit2=su/10; suit1=suite-10*suit2; /*suit2=suite/10; suit1=suite-10*suit2;*/ /*if (su3) fprintf(FTRACE,"DEBUT ISOLVE icout=%d\\\\\n",ICOUT);*/ if ((r=init_janus(II,XX,suit1))) return(r); stopc=0; if (su3) stopc=9; if (suit2) VRESULT=iprimal(XX,XX->minimum,ICOUT,stopc); /*stopcout*/ else VRESULT=is2(XX,&AVV,&ARR); /*ssss*/ if (VISU>=ZVL) vzlast(XX) ; return(VRESULT) ; }