pip.c

Go to the documentation of this file.

/*
 
  $Id: pip.c 23065 2016-03-02 09:05:50Z coelho $
 
  Copyright 1989-2016 MINES ParisTech
 
  This file is part of PIPS.
 
  PIPS is free software: you can redistribute it and/or modify it
  under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  any later version.
 
  PIPS is distributed in the hope that it will be useful, but WITHOUT ANY
  WARRANTY; without even the implied warranty of MERCHANTABILITY or
  FITNESS FOR A PARTICULAR PURPOSE.
 
  See the GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with PIPS.  If not, see <http://www.gnu.org/licenses/>.
 
*/
#ifdef HAVE_CONFIG_H
    #include "pips_config.h"
#endif
/* Name     : pip.c
 * Package  : pip
 * Author   : F. Dumontet, A. Platonoff and A. Leservot
 * Date     : 30 july 1993
 * Historic :
 * - 16 nov 93, few changes (e.g. "test1" to "pip_in", etc.), AP
 * - 08 dec 93, creation of a new package named pip.
 * - 10 dec 93, new version of pip_solve : direct call to pip. AL
 * - 31 jul 95, cleaning, pip_solve_min, pip_solve_max, old_pip_solve,
 * old2_pip_solve, AP
 *
 * Documents:
 * Comments : file containing the functions that resolve a Parametric Integer
 * Programming problem with the PIP solver.
 */
 
/* Ansi includes        */
#include <setjmp.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <stdlib.h>
 
/* Newgen includes      */
#include "genC.h"
 
/* C3 includes          */
#include "boolean.h"
#include "arithmetique.h"
#include "vecteur.h"
#include "contrainte.h"
#include "ray_dte.h"
#include "sommet.h"
#include "sg.h"
#include "sc.h"
#include "polyedre.h"
#include "matrix.h"
 
/* Pips includes        */
#include "ri.h"
#include "constants.h"
#include "ri-util.h"
#include "misc.h"
#include "bootstrap.h"
#include "complexity_ri.h"
#include "database.h"
#include "graph.h"
#include "dg.h"
#include "paf_ri.h"
#include "parser_private.h"
#include "property.h"
#include "reduction.h"
#include "text.h"
#include "paf-util.h"
#include "pip.h"
 
/* Macros and functions Used by the old_pip_solve version       */
#define PIP_BIN "pip4"
#define PIP_OPTION "-s"
#define PIP_IN_FILE "pip_in"
#define PIP_OUT_FILE "pip_out"
 
/* Variables for the direct call to PIP version */
#define INLENGTH 1024
long    int cross_product, limit;/* External variables for direct call to PIP */
int     allocation, comptage;
char    inbuff[INLENGTH];
int     inptr = 256;
int     proviso = 0;
int     verbose = 0; /* Should not be used : put here for Pip copatibility */
FILE    *dump = NULL;
char    dump_name[] = "XXXXXX";
 
/* Global variables     */
quast   quast_act;
Pbase   base_var_ref, base_ref, old_base, old_base_var;
int     ind_min_max;
 
/*===========================================================================*/
/* quast old_pip_solve(Psysteme ps_dep, ps_context, int nb_unknowns,
 *                 int min_or_max): Pip resolution.
 *
 * Parameters: ps_dep:     system of constraints on the unknowns
 *                         (and parameters)
 *             ps_context: system of constraints on the parameters
 *             nb_unknowns: number of unknowns in "ps_dep"
 *             min_or_max:  tag for min or max resolution
 *
 * Result: a quast giving the value of the unknowns
 *
 * Note: The basis of "ps_dep" must contain all the unknowns and the
 * parameters, with first the unknowns in order to allow us to catch them using
 * "nb_unknowns".
 */
/* Obsolete, AP July 31st 95 :
   
   quast old_pip_solve(ps_dep, ps_context, nb_unknowns, min_or_max)
   Psysteme ps_dep, ps_context;
   int nb_unknowns;
   int min_or_max;
   {
     Pvecteur pvect;
     int aux, n, i, res, infinite_num;
     char *com = "essai";
   
     if(min_or_max == PIP_SOLVE_MIN) {
                                       ind_min_max = 0;
                                       infinite_num = -1;
                                     }
     else {
            ps_dep = converti_psysmin_psysmax(ps_dep, nb_unknowns);
            ind_min_max = 1;
            infinite_num = vect_size((ps_dep)->base);
          }
   
     base_var_ref = base_dup((ps_dep)->base);
     old_base_var = base_dup(base_var_ref);
     for(i = 1, pvect = base_var_ref; i<nb_unknowns; i++) { pvect = pvect->succ; }
     base_ref = pvect->succ;
     pvect->succ = NULL;
     old_base = base_dup(base_ref);
   
     if(! SC_EMPTY_P(ps_context))
       ps_context->base = base_dup(base_ref);
   
     res = ecrit_probleme2(com, ps_dep, ps_context, nb_unknowns,
                           infinite_num, vect_size(base_ref));
 
     n = 1;
     while(n & 0xff) {
                       int m;
                       m = fork();
                       while(m == -1) {
                                        char answer;
   
                                        fprintf(stderr, "Fork failed for PIP\n\t Do you want to retry (y/n) ?\n");
                                        scanf("%c\n", &answer);
                                        if(answer == 'y')
                                          m = fork();
                                        else
                                          exit(1);
                                        fprintf(stderr, "\n");
                                      }
                       if(m == 0) {
                                    execlp(PIP_BIN, PIP_BIN, PIP_OPTION, PIP_IN_FILE, PIP_OUT_FILE, (char *) 0);
                                  }
 
                       wait(&n);
                     }
 
     if((quayyin = fopen(PIP_OUT_FILE,"r")) == NULL) {
                                                       fprintf(stderr, "Cannot open file %s\n", PIP_OUT_FILE);
                                                       exit(1);
                                                     }
     aux = quayyparse();
   
     return(quast_act);
   }
*/
 
/*===========================================================================*/
/* quast old2_pip_solve(Psysteme ps_dep, ps_context, int nb_unknowns,
 *                 int min_or_max): Pip resolution.  AL 6/12/93
 *
 * Parameters: ps_dep:     system of constraints on the unknowns
 *                         (and parameters)
 *             ps_context: system of constraints on the parameters
 *             nb_unknowns: number of unknowns in "ps_dep"
 *             int_or_rat:      integer or rational resolution.
 *             min_or_max:  tag for min or max resolution
 *
 * Result: a quast giving the value of the unknowns
 *
 * Note: The basis of "ps_dep" must contain all the unknowns and the
 * parameters, with first the unknowns in order to allow us to catch them using
 * "nb_unknowns".
 *
 * Call to PIP is direct : no input or output file are produced.
 */
/* Obsolete, AP July 31st 95 :
 
quast old2_pip_solve(ps_dep, ps_context, nb_unknowns, min_or_max)
Psysteme ps_dep, ps_context;
int nb_unknowns, min_or_max;
{
        Pvecteur pvect;
        int     i, infinite_num;
 
        int     q, bigparm, ni, nvar, nq, non_vide, nparm, nc, p, xq;
        char*   g;
        Tableau *ineq, *context, *ctxt;
 
 
        debug_on("PIP_DEBUG_LEVEL");
        debug(5, "pip_solve", "begin\n");
        if (get_debug_level()>4) {
                debug(5, "pip_solve", "Input Psysteme:\n");
                fprint_psysteme( stderr, ps_dep );
                debug(5, "pip_solve", "Input Context:\n");
                fprint_psysteme( stderr, ps_context );
                fprintf(stderr, "Number of variables : %d\n", nb_unknowns);
        }
 
        vect_erase_var(&(ps_dep->base), (Variable) TCST);
        ps_dep->dimension = vect_size((ps_dep)->base);
 
        if(min_or_max == PIP_SOLVE_MIN) {
                ind_min_max = 0;
                infinite_num = -1;
        }
        else {
                ps_dep = converti_psysmin_psysmax(ps_dep, nb_unknowns);
 
                vect_erase_var(&(ps_dep->base), (Variable) TCST);
                ps_dep->dimension = vect_size((ps_dep)->base);
 
                ind_min_max = 1;
                infinite_num = vect_size((ps_dep)->base);
        }
 
        base_var_ref = base_dup((ps_dep)->base);
        old_base_var = base_dup(base_var_ref); 
 
        for(i = 1, pvect = base_var_ref; i<nb_unknowns; i++)
                { pvect = pvect->succ; }
        base_ref = pvect->succ;                 
        pvect->succ = NULL;                     
        old_base = base_dup(base_ref);
 
        if(! SC_EMPTY_P(ps_context) ) {
                ps_context->base = base_dup(base_ref);
                ps_context->dimension = vect_size( base_ref );
        }
 
        nvar    = nb_unknowns;
        nparm   = ps_dep->dimension - nb_unknowns;
        ni      = (ps_dep->nb_eq * 2) + ps_dep->nb_ineq;;
        nc      = ((ps_context == NULL)? 0 :
                         (ps_context->nb_eq * 2) + ps_context->nb_ineq );
        bigparm = infinite_num;
        nq      = 1;
        debug(5, "pip_solve", "%d  %d  %d  %d  %d  %d\n",
                                nvar, nparm, ni, nc, bigparm, nq );
 
        limit   = 0L;
        sol_init();
        tab_init();
        cross_product = 0;
        g       = tab_hwm();
        ineq    = sc_to_tableau(ps_dep, nb_unknowns);
        if (ps_context != NULL) {
                context = sc_to_tableau(ps_context, 0);
        }
        else context = NULL;
        xq = p = sol_hwm();
 
        if (nc) {
                ctxt = expanser(context, nparm, nc, nparm+1, nparm, 0, 0);
                traiter( ctxt, NULL, True, UN, nparm, 0, nc, 0, -1 );
                non_vide = is_not_Nil(p);
                sol_reset(p);
        }
        else non_vide = True;
        if ( non_vide ) {
                traiter( ineq, context, nq, UN, nvar, nparm, ni, nc, bigparm );
                q = sol_hwm();
                init_new_base();
                while((xq = new_sol_edit(xq)) != q);
                sol_reset(p);
        }
        else quast_act = quast_undefined;
        tab_reset(g);
 
        if (get_debug_level()>5) {
                imprime_quast( stderr, quast_act );
        }
        debug_off();
        return(quast_act);
}
*/
 
 
/*=======================================================================*/
/* quast pip_solve_min_with_big(Psysteme ps_dep, ps_context, int nb_unknowns,
 *                 char *big): Pip resolution.     
 *
 * Parameters: ps_dep:     system of constraints on the unknowns
 *                         (and parameters)
 *             ps_context: system of constraints on the parameters
 *             nb_unknowns: number of unknowns in "ps_dep"
 *             big: big parameter 's name.
 *
 * Result: a quast giving the value of the unknowns
 *
 * Note: The basis of "ps_dep" must contain all the unknowns and the
 * parameters, with first the unknowns in order to allow us to catch them
 * using "nb_unknowns".
 *
 * Call to PIP is direct : no input or output file are produced.
 *
 * Note: this function is the same as pip_solve_min() but it has a big
 * parameter given by the use.
 */
 
quast pip_solve_min_with_big(ps_dep, ps_context, pv_unknowns, big)
 
 Psysteme   ps_dep, ps_context;
 Pvecteur   pv_unknowns;
 char       *big; 
{
 /* FD variables */
 int        infinite_num;
 
 /* Pip variables */
 int        q, bigparm, ni, nvar, nq, non_vide, nparm, nc, p, xq;
 char*      g;
 Tableau    *ineq, *context, *ctxt;
 
 /* AC variables for the big parameter */
 bool    not_found = true;
 list       lbase;
 entity     ent;
 int        nb_unknowns, len_big;
 
 debug_on("PIP_DEBUG_LEVEL");
 debug(5, "pip_solve_min_with_big", "begin\n");
 /* trace file */
 nb_unknowns = vect_size( pv_unknowns ); 
 if (get_debug_level()>4) 
    {
     debug(5, "pip_solve_min_with_big", "Input Psysteme:\n");
     fprint_psysteme( stderr, ps_dep );
     debug(5, "pip_solve_min_with_big", "Input Context:\n");
     fprint_psysteme( stderr, ps_context );
     fprintf(stderr, "Number of variables : %d\n", nb_unknowns);
    }
 
 /* We set the env for the kind of resolution desired (i.e. Min) */
 /* and we get the order of the bi parameter called "big" in the */
 /* base of ps_dep                                               */
 
 /* Normalize base of system ps_dep */
 pv_unknowns = base_normalize( pv_unknowns );
 vect_erase_var( &ps_dep->base, TCST );
 ps_dep->base = vect_add_first(pv_unknowns, ps_dep->base);
 ps_dep->dimension = vect_size( ps_dep->base );
 
 ind_min_max = 0;
 infinite_num = 1;
 len_big = strlen(big);
 lbase = base_to_list(ps_dep->base);
 while ((not_found) && (lbase != NIL))
    {
     ent = ENTITY(CAR(lbase));
     if (!strncmp(entity_local_name(ent), big, len_big))
         not_found = false; 
     else
        {
         infinite_num++;
         lbase = CDR(lbase);
        }
    }
 if (not_found) infinite_num = -1;
 
 /* Computation of the basis (unknowns and parameters)
  * base of ps_dep is the reference for all unknowns and parameters.
  */
 if (ps_context != NULL) {
     vect_erase_var( &ps_context->base, TCST);
     ps_context->dimension = vect_size( ps_context->base );
     base_var_ref = base_union( ps_dep->base, ps_context->base );
 }
 else base_var_ref = ps_dep->base;
 base_var_ref = vect_add_first(pv_unknowns, base_var_ref);
 base_ref = vect_substract(base_var_ref, pv_unknowns);
 
 old_base_var = base_dup(base_var_ref);  /* Total base of ps_dep */
 old_base = base_dup(base_ref);
 
 if (!SC_EMPTY_P(ps_context)) 
    {
     ps_context->base = base_dup(base_ref);
     ps_context->dimension = vect_size( base_ref );
    }
 
 /* Set PIP variables. Comes from ecrit_probleme2 */
 nvar  = nb_unknowns;
 nparm = ps_dep->dimension - nb_unknowns;
 ni    = (ps_dep->nb_eq * 2) + ps_dep->nb_ineq;;
 nc    = ((ps_context == NULL)? 0 :
                       (ps_context->nb_eq * 2) + ps_context->nb_ineq );
 bigparm = infinite_num;
 nq      = 1;
 debug(5, "pip_solve_min_with_big", "%d  %d  %d  %d  %d  %d\n",\
                        nvar, nparm, ni, nc, bigparm, nq );
 
 /* Prepare to call PIP */
 limit   = 0L;
 sol_init();
 tab_init();
 cross_product = 0;
 g    = tab_hwm();
 ineq = sc_to_tableau(ps_dep, nb_unknowns);
 if (ps_context != NULL) 
       context = sc_to_tableau(ps_context, 0);
 else  context = NULL;
 xq = p = sol_hwm();
 
 /* Verification de la non vacuite du contexte */
 if (nc)
    {
     ctxt = expanser(context, nparm, nc, nparm+1, nparm, 0, 0);
     traiter( ctxt, NULL, PIP_SOLVE_RATIONAL, UN, nparm, 0, nc, 0, -1 );
     non_vide = is_not_Nil(p);
     sol_reset(p);
    }
 else non_vide = True;
 
 if (non_vide) 
    {
     traiter(ineq, context, PIP_SOLVE_RATIONAL, UN, nvar, nparm, ni, nc, bigparm);
     q = sol_hwm();
     init_new_base();
     /* We read solution and put it in global quast_act */
     while((xq = rational_sol_edit(xq)) != q);
     sol_reset(p);
    }
 else quast_act = quast_undefined;
 tab_reset(g);
 
 if (get_debug_level()>5) imprime_quast(stderr, quast_act);
                
 debug_off();
 
 return(quast_act);
}
 
/*=======================================================================*/
/* Pvecteur vect_add_first( pv1, pv2 )                          AL 16 02 94
 * 
 * Suppress all variables of pv1 in pv2 and add pv1 before pv2.
 * Keep pv1 order, not pv2's order!
 */
Pvecteur vect_add_first( pv1, pv2 )
Pvecteur pv1, pv2;
{
        Pvecteur pv, pv11, pv22;
 
        if (pv1 == NULL) return pv2;
        pv11 = vect_reversal( vect_dup( pv1 ) );
        pv22 = vect_dup( pv2 );
 
        for( pv = pv11 ; pv != NULL; pv = pv->succ) {
                vect_erase_var( &pv22, pv->var );
                if (pv->succ != NULL) continue;
                pv->succ = pv22;
                break;
        }
        return pv11;
}
 
/*===========================================================================*/
/* quast pip_solve(Psysteme ps_dep, ps_context, Pvecteur pv_unknowns,
 *                 int int_or_rat, int min_or_max): Pip resolution.  AL 6/12/93
 *
 * Parameters: ps_dep:     system of constraints on the unknowns
 *                         (and parameters)
 *             ps_context: system of constraints on the parameters
 *             pv_unknowns: ordered vector of unknowns in "ps_dep"
 *             int_or_rat:      integer (1)  or rational (0)  resolution.
 *             min_or_max:  tag for min or max resolution
 *
 * Result: a quast giving the value of the unknowns
 *
 * Call to PIP is direct : no input or output file are produced.
 */
quast pip_solve(ps_dep, ps_context, pv_unknowns, int_or_rat, min_or_max)
Psysteme ps_dep, ps_context;
Pvecteur pv_unknowns;
int      int_or_rat, min_or_max;
{
        /* FD variables */
        int     infinite_num, nb_unknowns;
 
        /* Pip variables */
        int     q, bigparm, ni, nvar, nq, non_vide, nparm, nc, p, xq;
        char*   g;
        Tableau *ineq, *context, *ctxt;
 
 
        /* Initialization */
        debug_on("PIP_DEBUG_LEVEL");
        debug(5, "pip_solve", "begin\n");
        if (ps_dep == NULL) {
                user_warning("new_pip_solve", "\nInput Psysteme is empty !\n");
                return quast_undefined;
        }
        nb_unknowns = vect_size( pv_unknowns ); 
 
        /* trace file */
        if (get_debug_level()>4) {
                debug(5, "pip_solve", "Input Psysteme:\n");
                fprint_psysteme( stderr, ps_dep );
                debug(5, "pip_solve", "Input Context:\n");
                fprint_psysteme( stderr, ps_context );
                fprintf(stderr, "Number of variables : %d\n", nb_unknowns);
        }
 
 
        /* Normalize base of system ps_dep */
        pv_unknowns = base_normalize( pv_unknowns );
        vect_erase_var( &ps_dep->base, TCST );
        ps_dep->base = vect_add_first(pv_unknowns, ps_dep->base);
        ps_dep->dimension = vect_size( ps_dep->base );
 
 
        /* We set the env for the kind of resolution desired (Min or Max) */
        if(min_or_max == PIP_SOLVE_MIN) {
                ind_min_max = 0;
                infinite_num = -1;
        }
        else {
                ps_dep = converti_psysmin_psysmax(ps_dep, nb_unknowns);
                ind_min_max = 1;
                infinite_num = vect_size((ps_dep)->base);
        }
 
        /* Computation of the basis (unknowns and parameters)
         * base of ps_dep is the reference for all unknowns and parameters.
         */
        if (ps_context != NULL) {       
                vect_erase_var( &ps_context->base, TCST);
                ps_context->dimension = vect_size( ps_context->base );
                base_var_ref = base_union( ps_dep->base, ps_context->base );
        }
        else base_var_ref = ps_dep->base;
        base_var_ref = vect_add_first(pv_unknowns, base_var_ref);
        base_ref = vect_substract(base_var_ref, pv_unknowns);
 
        old_base_var = base_dup(base_var_ref);  /* Total base of ps_dep */
        old_base = base_dup(base_ref);
 
        if(! SC_EMPTY_P(ps_context) ) {
                ps_context->base = base_dup(base_ref);
                ps_context->dimension = vect_size( base_ref );
        }
 
        /* Set PIP variables. Comes from ecrit_probleme2 */
        nvar    = nb_unknowns;
        nparm   = ps_dep->dimension - nb_unknowns;
        ni      = (ps_dep->nb_eq * 2) + ps_dep->nb_ineq;;
        nc      = ((ps_context == NULL)? 0 :
                         (ps_context->nb_eq * 2) + ps_context->nb_ineq );
        bigparm = infinite_num;
        nq      = int_or_rat;
        debug(5, "pip_solve", "%d  %d  %d  %d  %d  %d\n",
                                nvar, nparm, ni, nc, bigparm, nq );
 
        /* Prepare to call PIP */
        limit   = 0L;
        sol_init();
        tab_init();
        cross_product = 0;
        g       = tab_hwm();
        ineq    = sc_to_tableau(ps_dep, nb_unknowns);
        if (ps_context != NULL) {
                context = sc_to_tableau(ps_context, 0);
        }
        else context = NULL;
        xq = p = sol_hwm();
 
        /* Verification de la non vacuite du contexte */
        if (nc) {
                ctxt = expanser(context, nparm, nc, nparm+1, nparm, 0, 0);
                traiter( ctxt, NULL, nq, UN, nparm, 0, nc, 0, -1 );
                non_vide = is_not_Nil(p);
                sol_reset(p);
        }
        else non_vide = True;
        if ( non_vide ) {
                traiter( ineq, context, nq, UN, nvar, nparm, ni, nc, bigparm );
                q = sol_hwm();
                init_new_base();
                /* We read solution and put it in global quast_act */
                if( int_or_rat == PIP_SOLVE_INTEGER ) 
                        while((xq = integer_sol_edit(xq)) != q);
                else while((xq = rational_sol_edit(xq)) != q);
                sol_reset(p);
        }
        else quast_act = quast_undefined;
        tab_reset(g);
 
        if (get_debug_level()>5) {
                imprime_quast( stderr, quast_act );
        }
        debug_off();
        return(quast_act);
}
        
 
/*===================================================================*/
quast pip_integer_min( ps_dep, ps_context, pv_unknowns )
Psysteme        ps_dep, ps_context;
Pvecteur        pv_unknowns;
{
        return( pip_solve( ps_dep, ps_context, pv_unknowns, 
                        PIP_SOLVE_INTEGER, PIP_SOLVE_MIN ));
}
 
/*===================================================================*/
quast pip_integer_max( ps_dep, ps_context, pv_unknowns )
Psysteme        ps_dep, ps_context;
Pvecteur        pv_unknowns;
{
        return( pip_solve( ps_dep, ps_context, pv_unknowns, 
                        PIP_SOLVE_INTEGER, PIP_SOLVE_MAX ));
}
 
/*===================================================================*/
quast pip_rational_min( ps_dep, ps_context, pv_unknowns )
Psysteme        ps_dep, ps_context;
Pvecteur        pv_unknowns;
{
        return( pip_solve( ps_dep, ps_context, pv_unknowns, 
                        PIP_SOLVE_RATIONAL, PIP_SOLVE_MIN ));
}
 
/*===================================================================*/
quast pip_rational_max( ps_dep, ps_context, pv_unknowns )
Psysteme        ps_dep, ps_context;
Pvecteur        pv_unknowns;
{
        return( pip_solve( ps_dep, ps_context, pv_unknowns, 
                        PIP_SOLVE_RATIONAL, PIP_SOLVE_MAX ));
}
 
/*===================================================================*/
/* void pip_solve_min(Psysteme ps_dep, ps_context, int nb_unknowns):
 * Pip resolution for the minimum.
 */
/* Obsolete, AP July 31st 95 :
   
   quast pip_solve_min(ps_dep, ps_context, nb_unknowns)
   Psysteme ps_dep, ps_context;
   int nb_unknowns;
   {
     return(old2_pip_solve(ps_dep, ps_context, nb_unknowns, PIP_SOLVE_MIN));
   }
   */
 
/*===================================================================*/
/* void pip_solve_max(Psysteme ps_dep, ps_context, int nb_unknowns):
 * Pip resolution for the maximum.
 */
/* Obsolete, AP July 31st 95 :
 
   quast pip_solve_max(ps_dep, ps_context, nb_unknowns)
   Psysteme ps_dep, ps_context;
   int nb_unknowns;
   {
   return(old2_pip_solve(ps_dep, ps_context, nb_unknowns, PIP_SOLVE_MAX));
   }
*/
 
/*===================================================================*/
 
/* Useful for the sorting of the variables in the system's Pvecteur.  
*/
static Pbase base_for_sort = (Pbase)NULL;
 
/* function for qsort. 
 * rank_of_variable returns 0 for TCST, hence the strange return
 */
static int compare_variables_in_base(pv1, pv2)
Pvecteur *pv1, *pv2;
{
    int 
        rank_1 = rank_of_variable(base_for_sort, var_of(*pv1)),
        rank_2 = rank_of_variable(base_for_sort, var_of(*pv2));
 
    message_assert("variable not in global vector", rank_1>=0 && rank_2>=0);
 
    /* TCST must be last, but given as 0
     */ 
    return((rank_1==0 || rank_2==0) ? rank_1-rank_2 : rank_2-rank_1);
}
 
/*===================================================================*/
/* void sort_psysteme(Psysteme ps, Pvecteur pv): sorts the system "ps"
 * according to the vector "pv". In order to avoid any side effects,
 * this vector is duplicated into a base "new_base". When the sorting
 * is done, this new base becomes the basis of this system.
 *
 * "pv" MUST contain all the variables that appear in "ps".
 * "pv" DOES NOT HAVE TO contain TCST.
 *
 * modified by FC, 29/12/94.
 * I removed the vect_tri call and the associated memory leaks.
 * sc_vect_sort called instead.
 * maybe the sc base is lost? It depends whether it is pv or not,
 * and what is done with pv by the caller. I would suggest a base_rm.
 */
void sort_psysteme(ps, pv)
Psysteme ps;
Pvecteur pv;
{
    if(SC_EMPTY_P(ps))  return;
 
    base_for_sort = base_dup(pv);
    sc_vect_sort(ps, compare_variables_in_base);
 
    ps->dimension = base_dimension(base_for_sort);
    ps->base = base_for_sort;
    base_for_sort = (Pbase)NULL;
}
 
/*************************************************************************/