Logo Search packages:      
Sourcecode: gambit version File versions  Download package

pelutils.cc

/*
  This file contains the implementation information that had been in globals.h,
and also the various implementation files that formerly resided in the Utils
subdirectory of the Pelican distribution.  
*/

#include "pelutils.h"


/**************************************************************************/
/*********************** implementations from Mem.c ***********************/
/**************************************************************************/

/*
   **  Mem.c               Created 1-27-95                Birk Huber
   **
   **  Memory Management for Pelican Shell. 
   **  
   **  Maintains storage for the basic unit of storage from which
   **  lists are built.
   ** 
   **  A node consists of left,and right (values,type) combinations, and
   **  a boolean for garbage collection. The values are often pointers to
   **  other structures wich must be freed when nodes are freed.
   **
   ** The Free_Stack is used to store a list of (known) free nodes, 
   ** The Local_Stack points to a list of pointers to nodes
   **                           which are known to be in use.
   **               
   **  Any time a new_node is requested and the Free_Stack is 
   **  empty garbage collection is initiated. Garbage collection 
   **  proceeds by  first marking all nodes accessable from the 
   **  Local_Stack and then putting the rest on the free list.
   **
   **  Any time garbage collection fails node_more_store() is called
   **  to add another block of nodes to the pool of available nodes.
   **
 */

#define BLOCKSIZE 10000
#define YES 0
#define NO 1

#define Node_LT(n) ((n->LT))
#define Node_L(n) ((n->L))
#define Node_RT(n) ((n->RT))
#define Node_R(n) ((n->R))
#define Node_Marked(n) ((n)->Mark)
#define Node_Unset_Mark(n)((n)->Mark=NO)
#define Node_Set_Mark(n)((n)->Mark=YES)

static int node_gc();
/* static int mark(); */
static int node_more_store();
static node node_free(node N);




struct node_t {
    int LT, RT, Mark;
    union {
      void *ptr;
      int ival;
      double dval;
    } L, R;
};

/* 
   ** 
 */
typedef struct node_block_t *node_block;
struct node_block_t {
    struct node_t store[BLOCKSIZE];
    struct node_block_t *next;
};



static int N_MALLOC = 0;      /*number of mallocs called from module */
static int N_FREE = 0;        /*number of frees called from this module */
static node_block Node_Store = 0;   /*beginning of node storage */
static node Free_List = 0;    /*top of free node list */
static local_v Locals_Stack = 0;    /*stack of  declared pointers 
                                 into node storage */



node node_set_ptr(node N, void *v, int tp, int side)
{
    if (N == 0)
      bad_error("null node to node_set_ptr");
    if (side == LEFT) {
      N->L.ptr = v;
      N->LT = tp;
    } else if (side == RIGHT) {
      N->R.ptr = v;
      N->RT = tp;
    } else
      bad_error("bad side given in node_set_ptr");
    return N;
}

node node_set_int(node N, int v, int tp, int side)
{
    if (N == 0)
      bad_error("null node to node_set_int");
    if (side == LEFT) {
      N->L.ival = v;
      N->LT = tp;
    } else if (side == RIGHT) {
      N->R.ival = v;
      N->RT = tp;
    } else
      bad_error("bad side given in node_set_int");
    return N;
}


node node_set_double(node N, double v, int tp, int side)
{
    if (N == 0)
      bad_error("null node to node_set_double");
    if (side == LEFT) {
      N->L.dval = v;
      N->LT = tp;
    } else if (side == RIGHT) {
      N->R.dval = v;
      N->RT = tp;
    } else
      bad_error("bad side given in node_set_double");
    return N;
}



int node_get_int(node N, int side)
{
    if (N == 0)
      bad_error("null node to node_get_int");
    if (side == LEFT)
      return N->L.ival;
    else if (side == RIGHT)
      return N->R.ival;
    bad_error("bad side given in node_get_int");
    return 0;
}
double node_get_double(node N, int side)
{
    if (N == 0)
      bad_error("null node to node_get_ptr");
    if (side == LEFT)
      return N->L.dval;
    else if (side == RIGHT)
      return N->R.dval;
    bad_error("bad side given in node_get_double");
    return 0;
}

void *node_get_ptr(node N, int side)
{
    if (N == 0)
      bad_error("null node to node_get_ptr");
    if (side == LEFT)
      return N->L.ptr;
    else if (side == RIGHT)
      return N->R.ptr;
    bad_error("bad side given in node_get_ptr");
    return 0;
}

int node_set_type(node N, int Tp, int side)
{
    if (N == 0)
      bad_error("null node to node_get_type");
    if (side == LEFT)
      return N->LT;
    else if (side == RIGHT)
      return N->RT=Tp;
    bad_error("bad side given in node_get_type");
    return 0;
}
int node_get_type(node N, int side)
{
    if (N == 0)
      bad_error("null node to node_get_type");
    if (side == LEFT)
      return N->LT;
    else if (side == RIGHT)
      return N->RT;
    bad_error("bad side given in node_get_type");
    return 0;
}


node Cons(node N1, node N2)
{
    node Res;
    LOCS(2);
    PUSH_LOC(N1);
    PUSH_LOC(N2);
    Res = node_new();
    node_set_ptr(Res, (void *) N1, NODE, LEFT);
    node_set_ptr(Res, (void *) N2, NODE, RIGHT);
    POP_LOCS();
    return Res;
}

node Car(node N1)
{
    if (N1 == 0)
        return 0;
    return (node) node_get_ptr(N1, LEFT);
}

node Cdr(node N1)
{
    if (N1 == 0)
        return 0;
    return (node) node_get_ptr(N1, RIGHT);
}

int node_atomp(node N1)
{
    
    if ((N1!=0) && (node_get_type(N1, LEFT)) == NODE)
        return FALSE;
    else
        return TRUE;
}

int node_nullp(node N1)
{
    if (N1 == 0)
        return TRUE;
    else
        return FALSE;
}






char *mem_strdup(char *str){
  /*   char *strdup(char *); CAN'T DECLARE BUILTINS UNDER C++ */
  N_MALLOC++;
  return strdup(str);
}
void *mem_malloc(int sz)
{
    N_MALLOC++;
    return (void *) malloc(sz);
}
void mem_free(void *ptr)
{
    N_FREE++;
    free(ptr);
}

static node node_free(node N)
{
    Node_LT(N) = NODE;
    Node_L(N).ptr = 0;
    Node_RT(N) = NODE;
    Node_R(N).ptr = Free_List;
    Free_List = N;
    return N;
}

/* 
   **  storage is stored in blocks of size BLOCKSIZEg. whenever
   **  node_new is called and garbage collection fails to produce 
   **  new nodes (because there is no more free storage) 
   **  node_more_store() is called. 
   **  It is responcible for 
   **     a) claiming another block of nodes with malloc
   **     b) adding this block to the storage list
   **     c) putting all the new blocks on the free list.
   **  finally if this is the first block to be created the locals stack
   **  must be initialized.
 */
static int node_more_store()
{
    int i;
    node_block next_node = 0;
  
    next_node = (node_block) malloc(sizeof(struct node_block_t));
    if (next_node == 0) {
#ifdef LOG_PRINT
      fprintf(stderr, "malloc failing in node_more_store\n")
#endif
;
      return FALSE;
    }
    next_node->next = Node_Store;
    Node_Store = next_node;
    for (i = 0; i < BLOCKSIZE; i++)
      Node_Unset_Mark(node_free(Node_Store->store + i));

    return TRUE;
}


int node_init_store()
{
    if (node_more_store() == FALSE)
      bad_error("node_init_store failing");
    return TRUE;
}



void node_free_store()
{
    int i;
    node_block junk;
    printf("in free_store()\n");

    while (Node_Store != 0) {
      /* free any space held by elements in the block */
      for (i = 0; i < BLOCKSIZE; i++)
          atom_free(Node_Store->store + i);
      Node_Store = (junk = Node_Store)->next;
      free((char *) junk);
    }
    printf("N_malloc = %d,  N_free=%d;=\n \n", N_MALLOC, N_FREE);
}


/* 
   ** Allocate nodes from free list-- If no nodes available initiate
   ** a garbage collection.  Note: routines which call new-node directly
   ** must protect their arguments, and localvariables themselve rather
   ** than make any assumptions about wheather they have been saved.
 */

node node_new()
{
    node res;
    if (Free_List != 0) {
      res = Free_List;
      Free_List = (node) node_get_ptr(Free_List, RIGHT);
      node_set_ptr(res, (void *) 0, NODE, LEFT);
      node_set_ptr(res, (void *) 0, NODE, RIGHT);
      return res;
    }
    if (node_gc() != 0)
      return node_new();
    if (node_more_store() != FALSE)
      return node_new();
#ifdef LOPG_PRINT  
  fprintf(stderr, "out of nodes: system not giving more memory\n")
#endif
;
    abort();
    return node_new();
}

/*
   ** void node_push_local(node *local)
   **   puts the address of a local variable onto the Locals_Stack stack 
   **   ( of starting points for garbage collection). 
   **
   ** void node_pop_local()
   **     removes one ptr off the Locals_Stack stack.
 */

void print_locals_stack(void)
{
    local_v ptr;
    ptr = Locals_Stack;
    printf("printing local stack \n");
    while (ptr != 0) {
      node_print(*(ptr->val));
      ptr = ptr->next;
    }
    printf(" done \n");
}

void node_push_local(local_v loc, node * val)
{
#ifdef MEM_DEBUG
  if (Locals_Stack == loc)
    bad_error("trying to make circular stack\n");
#endif

    loc->next = Locals_Stack;
    loc->val = val;
    Locals_Stack = loc;
}

void node_pop_local()
{
    if (Locals_Stack == 0)
      bad_error("poping empty local stack");
    Locals_Stack = Locals_Stack->next;
}




/*
   **  Mark all nodes accessable from a node n.
   **  if n is marked Yes then both its subtrees are already 
   **  marked and should be skipped. 
   **  otherwise cal mark on any pointers to nodes stored in the 
   **  node (is_node is used to check this)
 */
static int mark(node n)
{
    if ((n != 0) && (Node_Marked(n) == NO)) {
      Node_Set_Mark(n);
      if (Node_LT(n) == NODE && Node_L(n).ptr != 0)
          mark((node) Node_L(n).ptr);
      else if (Node_LT(n) == NPTR && Node_L(n).ptr != 0)
          mark(*((node *) Node_L(n).ptr));
      if (Node_RT(n) == NODE && Node_R(n).ptr != 0)
          mark((node) Node_R(n).ptr);
      else if (Node_RT(n) == NPTR && Node_R(n).ptr != 0)
          mark(*((node *) Node_R(n).ptr));
      return 0;
    }
    return 1;
}

/*
   ** Collect garbage
 */
static int node_gc()
{
    int i, free_cnt = 0;
    node_block next_block = Node_Store;
    local_v ptr;
#ifdef LOG_PRINT   
  fprintf(stdout, "Colecting Garbage .. ")
#endif
;
    fflush(stdout);

    /*call mark on variables in local list */
    ptr = Locals_Stack;
    while (ptr != 0) {
      mark(*(ptr->val));
      ptr = ptr->next;
    }

    /* now put anything that is not already marked on the free_list */
    while (next_block != 0) {
      for (i = 0; i < BLOCKSIZE; i++) {
          if (Node_Marked(next_block->store + i) == YES)
            Node_Unset_Mark(next_block->store + i);
          else {
            atom_free(next_block->store + i);
            node_free(next_block->store + i);
            free_cnt++;
          }
      }
      next_block = next_block->next;
    }

#ifdef LOG_PRINT   
 fprintf(stdout, "Done  %d nodes freed\n", free_cnt)
#endif
;
    fflush(stdout);
    return free_cnt;
}


/**************************************************************************/
/********************** implementations from error.c **********************/
/**************************************************************************/

#ifdef GAMBIT_EXCEPTIONS
ErrorInPelican::~ErrorInPelican() { }

std::string ErrorInPelican::GetDescription(void) const
{
  return "Error somewhere in Pelican";
}
#endif

void bad_error(char *m)     /* generates an error message and aborts*/
{
#ifdef GAMBIT_EXCEPTIONS
  throw ErrorInPelican();
#endif

 #ifdef LOG_PRINT
       fprintf(stderr /* was Pel_Err */,"%s\n",m)
#endif
;
        exit(1);
}

void warning(char *m){
#ifdef LOG_PRINT
fprintf(stderr /* was Pel_Err */,"Warning:%s\n",m)
#endif
;}


/**************************************************************************/
/*********************** implementations from Rand.c **********************/
/**************************************************************************/

void srand48(long int seedval);

// WARNING: RDM added the following just to get to compile under BCC
//            I have no idea if this is correct!!!  
#if !defined(HAVE_SRAND48)
void srand48(long int seedval)
{
  srand(seedval);
}
#endif // HAVE_SRAND48
/*
** rand_seed  -- seed the random number generator with seedval.
*/
void rand_seed(long int seedval){ srand48(seedval);}

/*
**rand_int  -- return a random integer r with low<=r<=high
**             it produces double between low and high and rounds
**              by adding .5-epsilon and truncating, (if we just add
**              .5 then there is a slight chance we could end up with
**              high+1.)
*/
int rand_int(int low, int high){
      return (int)(low+drand48()*(high-low)+.499999999999); 
}



/*
**rand_double  -- return a random integer r with low<=r<=high
*/
double rand_double(int low, int high){
      return (drand48()*(high-low)+low);
}


/**************************************************************************/
/********************** implementations from Dlist.c **********************/
/**************************************************************************/

void Dlist_empty(node L){ Dlist_first(L)=0;}

/*
** invariants: Dnode_next(Dnode_prev(pos)) ==pos 
**             for all pointers to Dlist node entrees.
**             (NOTE: Dnode_next(L) := Dlist_first(L), 
**                    where L is list header) 
**
**             Dnode_prev(Dnode_next(pos))=pos
**             for all non-zero pointers to Dlist node entrees 
**             and also for list header
**             
*/

node Dlist_add(node L, node data){
  node tmp=0;
  LOCS(2);
  PUSH_LOC(data);
  PUSH_LOC(tmp);
  tmp=node_new();
  Dnode_link(tmp)=node_new();
  Dnode_data(tmp)=data;
  Dnode_prev(tmp)=L;
  Dnode_next(tmp)=Dlist_first(L);
  if (Dnode_next(tmp)!=0) Dnode_prev(Dnode_next(tmp))=tmp;
  Dlist_first(L)=tmp; 
  POP_LOCS();
  return tmp;
}


node Dlist_del(node L, node pos){
  if (Dnode_next(pos)!=0) Dnode_prev(Dnode_next(pos))=Dnode_prev(pos);
  if (Dnode_prev(pos)!=L) Dnode_next(Dnode_prev(pos))=Dnode_next(pos);
  else Dlist_first(L)=Dnode_next(pos);
  return (node)Dnode_data(pos);
}

node Dlist_data(node pos){ return (node)Dnode_data(pos);}

/* end Dlist.c */


/**************************************************************************/
/******************** implementations from Dmatrix.c **********************/
/**************************************************************************/

/*--------------------------------------------------------------- 
  vector/matrix type  a linear array of int, whith auxilary info.
       *) the number of elements that can be stored is in elt[0]
       *) the current number of rows is in elt[1]
       *) the current number of collumbs is in elt[2]
The actual data are then stored in row major order from elt[3] on
---------------------------------------------------------------*/
#ifndef DMATRIX_FAST
struct Dmatrix_t {
    int store;
    int nrows;
    int ncols;
    double *coords;
};
#endif

/*-------------------------------------------------------------
 vector access macroes (which ignore any rows except for first)
-------------------------------------------------------------*/
#define Vstore(V)  ((V->store))     /* maximum #elts available */
#define Vlength(V) ((V->nrows))     /* actual #elts stored      */
#define Vref1(V,i) (((V->coords)[i-1]))      /*acces ith elt (starting at 1) */
#define Vref0(V,i)  (((V->coords)[i]))    /*acces ith elt (starting at 0) */
#define Vref(V,i)  Vref1(V,i)

/*------------------------------------------------------------
 matrix access macroes
-------------------------------------------------------------*/
#define Mstore(V)  ((V->store))     /* maximum #elts available */
#define MMrows(V)  ((V->store/V->ncols))  /* maximum #rows          */
#define Mrows(V) ((V->nrows)) /* number rows stored */
#define Mcols(V) ((V->ncols)) /* number cols stored */
#define Mref1(V,i,j)(((V->coords)[(i-1)*(V->ncols)+j-1]))
#define Mref0(V,i,j)(((V->coords)[i*(V->ncols)+j]))
#define Mref(V,i,j)  Mref1((V),i,j)

#ifndef DMATRIX_FAST
double DMstore(Dmatrix M)
{
    return Mstore(M);
}
double DMMrows(Dmatrix M)
{
    return MMrows(M);
}
double DMrows(Dmatrix M)
{
    return Mrows(M);
}
double DMcols(Dmatrix M)
{
    return Mcols(M);
}
double *DMelts(Dmatrix M)
{
    return M->coords;
}
double *DMref_P(Dmatrix M, int i, int j)
{
    return &(Mref0(M, i, j));
}

#endif 

/*
   **   Constructor/Destructors for Dmatrixes
   ** 
   ** Dmatrix Dmatrix_free(int r, int c); 
   **       New Dmatrix cabable of holding r rows, and c collumbs.
   ** Dmatrix Dmatrix_new(Dmatrix V);
 */

Dmatrix Dmatrix_new(int r, int c)
{
    Dmatrix V;
    /*    void *mem_malloc(int); */
    V = (Dmatrix) mem_malloc(sizeof(struct Dmatrix_t));
    if (!V)
      bad_error("allocation failure in Dmatrix_new()");
    V->coords = (double *) mem_malloc(r * c * sizeof(double));
    if (!V)
      bad_error("allocation failure 2 in Dmatrix_new()");
    Mstore(V) = r * c;
    Mrows(V) = r;
    Mcols(V) = c;
    return V;
}

void Dmatrix_free(Dmatrix V)
{   /* void mem_free(); */
    if (V != 0 && V->coords != 0)
      mem_free((char *) (V->coords));
    if (V != 0)
      mem_free((char *) (V));
}
#undef mem_malloc
#undef mem_free
/*
   ** Dmatrix_resize(R,r,c)
   **   if R has enough storage to hold an rxc matrix resets
   **   row and columb entrees of r to r and c. otherwise
   **   frees R and reallocates an rxc matrix
 */
Dmatrix Dmatrix_resize(Dmatrix R, int r, int c)
{

    if (R == 0 || Mstore(R) < (r * c)) {
      if (R != 0)
          Dmatrix_free(R);
      R = Dmatrix_new(r, c);
    } else {
      Mrows(R) = r;
      Mcols(R) = c;
    }
    return R;
}


/*
   **  Dmatrix_print(M):  print an Dmatrix
   **    if M is null print <<>> and return fail.
   **    otherwise print matrix and return true.
 */
Dmatrix Dmatrix_fprint(FILE *fout, Dmatrix M)
{
    int i, j;

    if (M == 0) {
#ifdef LOG_PRINT
      fprintf(fout,"<<>>")
#endif
;
      return 0;
    }
#ifdef LOG_PRINT
    fprintf(fout,"<")
#endif
;
    for (i = 1; i <= Mrows(M); i++) {
#ifdef LOG_PRINT
      fprintf(fout,"<")
#endif
;
      for (j = 1; j <= Mcols(M); j++) {
#ifdef LOG_PRINT
          fprintf(fout,"%8.4f", Mref(M, i, j))
#endif
;
          if (j < Mcols(M))
#ifdef LOG_PRINT  
      fprintf(fout,", ")
#endif
;
      }
#ifdef LOG_PRINT
      fprintf(fout,">")
#endif
;
      if (i < Mrows(M))
#ifdef LOG_PRINT
          fprintf(fout,"\n")
#endif
;
    }
#ifdef LOG_PRINT
    fprintf(fout,">")
#endif
;
    return M;
}

/*
   ** matrix_add(M1,M2,&M3) -- Add two Dmatrixes:
   **  if M1, and M2 are incompatable (or null) complain and return false.
   **  if *M3 has too little storage (or is null) free *M3 if nescesary
   **                                            and create new storage.
 */
Dmatrix Dmatrix_add(Dmatrix M1, Dmatrix M2, Dmatrix * M3)
{
    int i, j;
    Dmatrix R = *M3;

    if (M1 == 0 || M2 == 0 || Mrows(M1) != Mrows(M2) || Mcols(M1) != Mcols(M2)) {
#ifdef LOG_PRINT
      fprintf(stderr, "matrix_add: dimensions dont match\n")
#endif
;
      return 0;
    }
    Dmatrix_resize(R, Mrows(M1), Mcols(M1));
    for (i = 1; i <= Mrows(M1); i++)
      for (j = 1; j <= Mcols(M1); j++)
          Mref(R, i, j) = Mref(M1, i, j) + Mref(M2, i, j);
    M3 = &R;
    return R;
}


/*
   **Dmatrix_mull(M1,M2,&M3) -- multiply two Dmatrixes:
   ** if M1, and M2 are incompatable (or null) complain and return false.
   ** if *M3 has too little storage (or is null) free *M3 if nescesary
   **                                           and create new storage.
   ** NOT USED YET
 */

Dmatrix Dmatrix_mull(Dmatrix M1, Dmatrix M2, Dmatrix * M3)
{
    int i, j, k;
    Dmatrix R = *M3;

    if (M1 == 0 || M2 == 0 || Mcols(M1) != Mrows(M2)) {
#ifdef LOG_PRINT
      fprintf(stderr, "Dmatrix_mull: incompatible matrices\n")
#endif
;
      return 0;
    }
    Dmatrix_resize(R, Mrows(M1), Mcols(M2));
    for (i = 1; i <= Mrows(M1); i++)
      for (j = 1; j <= Mcols(M2); j++) {
          Mref(R, i, j) = 0;
          for (k = 1; k <= Mcols(M1); k++)
            Mref(R, i, j) += Mref(M1, i, k) * Mref(M2, k, j);
      }
    M3 = &R;
    return R;
}

/*
   ** Dmatrix_dot(M1,M2,&M3) -- calculate M1*Transpose(M2):                
   ** if M1, and M2 are incompatable (or null) complain and return false.
   ** if *M3 has too little storage (or is null) free *M3 if nescesary 
   **                                           and create new storage.   
 */
Dmatrix Dmatrix_dot(Dmatrix M1, Dmatrix M2, Dmatrix M3)
{
    int i, j, k;

    if (M1 == 0 || M2 == 0 || Mcols(M1) != Mcols(M2)) {
#ifdef LOG_PRINT
      fprintf(stderr, "Dmatrix_dot: incompatible matrices\n")
#endif
;
      return 0;
    }
    M3 = Dmatrix_resize(M3, Mrows(M1), Mrows(M2));

    for (i = 1; i <= Mrows(M1); i++)
      for (j = 1; j <= Mrows(M2); j++) {
          Mref(M3, i, j) = 0;
          for (k = 1; k <= Mcols(M1); k++)
            Mref(M3, i, j) += Mref(M1, i, k) * Mref(M2, j, k);
      }
    return M3;
}


/*
** Dmatrix_equal(M1,M2)  
**     return true if matrices represented by M1 and M2 are equal.
**     false otherwise.                     NOT TESTED YET
*/
int Dmatrix_equal(Dmatrix M1, Dmatrix M2)
{
    int i, j;
    if (M1 == 0 && M2 == 0)
      return 1;
    if (M1==0||M2==0||(Mrows(M1)!=Mrows(M2))||(Mcols(M1)!=Mcols(M2)))
        return 0;
    for (i = 1; i <= Mrows(M1); i++)
      for (j = 1; j <= Mcols(M1); j++)
          if (Mref(M1, i, j) != Mref(M2, i, j))
            return 0;

    return 1;
}

/*
** Dmatrix_GQR
**    Use givens QR on a sqaure matrix:
**        after completion Q should be orthogonal, R triangular
**        and Q*A = original matrix A
*/
void Dmatrix_GQR(Dmatrix Q,Dmatrix A)
        /* A->ch<=A->rh=Q->rh=Q->ch */
   { int i,j,k;
     int r,c;
     double s,s1,s2;
     double t1,t2;
     r=Mrows(A);
     c=Mcols(A);
       for(i=1;i<=r;i++){
              for(k=1;k<=r;k++)Mref(Q,i,k)=0.0;
              Mref(Q,i,i)=1.0;}

       for (i=1;i<=c;i++)
         for (k=i+1;k<=r;k++)
              /* performing givens rotations to zero A[k][i] */
             if (Mref(A,k,i)!=0){
               s=sqrt(Mref(A,i,i)*Mref(A,i,i)+
                      Mref(A,k,i)*Mref(A,k,i));
               s1=Mref(A,i,i)/s;
               s2=Mref(A,k,i)/s;
               for(j=1;j<=c;j++) {
                          t1=Mref(A,i,j);
                          t2=Mref(A,k,j);
                          Mref(A,i,j)=s1*t1+s2*t2;
                          Mref(A,k,j)=-s2*t1+s1*t2;}
            /* actually doing givens row rotations on transpose Q */
               for(j=1;j<=r;j++){
                    t1=Mref(Q,j,i);
                    t2=Mref(Q,j,k);
                    Mref(Q,j,i)=s1*t1+s2*t2;
                    Mref(Q,j,k)=-s2*t1+s1*t2;}
             }
         }

/* 
** Dmatrix_Solve, Solve a square matrix in place
*/
#define LHS(i,j) Mref(LHS,i,j)
#define RHS(i) Vref(RHS,i)
void Dmatrix_Solve(Dmatrix LHS, Dmatrix RHS,int n){
 int i,j,k,m=n;
 double s,s1,s2,t1,t2;

 /*
 ** Use Givens rotations to triangularize LHS,i.e.LHS becomes Q*LHS.
 ** (triangular)
 ** and also apply same givens rotations to RHS i.e. RHS becomes Q*RHS.
 */
 for (i=1;i<=n;i++){
   for (k=i+1;k<=m;k++){
     if (LHS(k,i)!=0.0){
         /* compute rotation */
         s=sqrt(LHS(i,i)*LHS(i,i)+LHS(k,i)*LHS(k,i));
         s1=LHS(i,i)/s;
         s2=LHS(k,i)/s;
         /* apply rotation to LHS */
         for(j=1;j<=n;j++) {
             t1=LHS(i,j);
             t2=LHS(k,j);
             LHS(i,j)=s1*t1+s2*t2;
             LHS(k,j)=-s2*t1+s1*t2;
         }
         /* apply rotation to RHS */
         t1=RHS(i);
         t2=RHS(k);
         RHS(i)=s1*t1+s2*t2;
         RHS(k)=-s2*t1+s1*t2;
     }
   }
 }
 /*
 ** Back solve   R*X=Y
 **   R=top square portion of LHS. (2nx2nupper triangular)
 **   Y=top 2n entrees of RHS.
 **   results in X which solves original least squares problem.
 */
 for(j=n;j>=1;j--){
   t1=RHS(j);
   for(i=j+1;i<=n;i++) t1-=LHS(j,i)*RHS(i);
   t1=t1/LHS(j,j);
   RHS(j)=t1;
 }
}
#undef LHS
#undef RHS

/* end Dmatrix.c */

/**************************************************************************/
/********************* implementations from Imatrix.c *********************/
/**************************************************************************/

#ifndef min
#define min(i,j) ((i) < (j) ? (i): (j))
#endif
#ifndef FALSE
#define FALSE 0
#endif
#ifndef TRUE
#define TRUE 1
#endif


#ifndef IMATRIX_FAST 
/*--------------------------------------------------------------- 
  vector/matrix type  a linear array of int, whith auxilary info.
       *) the number of elements that can be stored is in elt[0]
       *) the current number of rows is in elt[1]
       *) the current number of collumbs is in elt[2]
The actual data are then stored in row major order from elt[3] on
---------------------------------------------------------------*/
struct Imatrix_t {
    int store;    /* maximum number of helts reserved */
    int topc;     /* effective number of columbs */
    int topr;     /* effective number of rows */
    int ncols;    /* number of elts between first elts of rows */
    int *elts;    /* ptr to first elt of block of space */
};

/*
**  a couple of functions which are defined in files that
**  I don't want to include right now, when things stableize
**      bad_error -- used mainly to stop things when a function
**                   is passed bad input, as things stablizes
**                   less catastrophic reactions may be introduced.
**      mem_alloc,
**      mem_free -- call malloc and free with some extra book-keeping.
*/
#endif

/*
void bad_error(char *);
char *mem_malloc(size_t);
void mem_free(char *);
*/

/*-------------------------------------------------------------
 vector access macroes (which ignore any rows except for first)
-------------------------------------------------------------*/

#define ImatrixVstore(V)  (((V)->store))  /* maximum #elts available */
#define ImatrixVlength(V) (((V)->topc))   /* actual #elts stored      */
#define ImatrixVref1(V,i)  (((V)->elts[i-1]))/* acces ith elt (starting at 1)*/
#define ImatrixVref0(V,i)  (((V)->elts[i]))/* acces ith elt (starting at 0) */
#define ImatrixVref(V,i)  ImatrixVref1(V,i)

/*------------------------------------------------------------
 matrix access macroes
-------------------------------------------------------------*/
#define ImatrixMstore(V)  (((V)->store))  /* maximum #elts available */
#define ImatrixMMrows(V)  (((V)->store/(V)->ncols))   /* maximum #rows    */
#define ImatrixMrows(V) (((V)->topr))     /* number rows stored */
#define ImatrixMcols(V) (((V)->topc))     /* number cols stored */
#define ImatrixMNcols(V) (((V)->ncols))   /* number cols stored */
#define ImatrixMref1(V,i,j) ((((V))->elts[(i-1)*((V)->ncols)+j-1]))
#define ImatrixMref0(V,i,j) (((V)->elts[(i*(V)->ncols)+j]))
#define ImatrixMref(V,i,j)  ImatrixMref1(V,i,j)


#ifndef IMATRIX_FAST
int IMstore(Imatrix M)
{
    return ImatrixMstore(M);
}
int IMMrows(Imatrix M)
{
    return ImatrixMMrows(M);
}
int IMrows(Imatrix M)
{
    return ImatrixMrows(M);
}
int IMcols(Imatrix M)
{
    return ImatrixMcols(M);
}
int *IMref1(Imatrix M, int i, int j)
{
    return &(ImatrixMref1(M, i, j));
}
#endif



/*
   **   Constructor/Destructors for Imatrixes
   ** 
   ** Imatrix Imatrix_free(int r, int c); 
   **       New Imatrix cabable of holding r rows, and c collumbs.
   ** Imatrix Imatrix_new(Imatrix V);
 */

Imatrix Imatrix_new(int r, int c)
{
    Imatrix V;
    V = (Imatrix) mem_malloc(sizeof(struct Imatrix_t));
    if (!V)
      bad_error("allocation failure 1 in Imatrix_new()");
    V->elts = (int *) mem_malloc(r * c * sizeof(int));
    if (!V)
      bad_error("allocation failure 2 in Imatrix_new()");
    V->store = r * c;
    V->topc = V->ncols = c;
    V->topr = r;
    return V;
}

void Imatrix_free(Imatrix V)
{
    if (V != 0) {
      if (V->elts != 0)
          mem_free((char *) (V->elts));
      mem_free((char *) (V));
    }
}

/*
   ** Imatrix_resize(R,r,c)
   **   Reset R to hold an r,by  c matrix.
   **   if R has enough storage to hold an rxc matrix resets
   **   row and columb entrees of r to r and c. otherwise
   **   frees R and reallocates an rxc matrix
   ** DOES NOT PRESERVE INDECIES OF EXISTING DATA
 */
Imatrix Imatrix_resize(Imatrix R, int r, int c)
{

    if (R == 0 || ImatrixMstore(R) < (r * c)) {
      if (R != 0)
          Imatrix_free(R);
      R = Imatrix_new(r, c);
    } else {
      ImatrixMrows(R) = r;
      ImatrixMcols(R) = c;
      ImatrixMNcols(R) = c;
    }
    return R;
}


Imatrix Imatrix_submat(Imatrix R, int r, int c)
{

    if (R == 0 || c > ImatrixMNcols(R) || r > ImatrixMrows(R) * ImatrixMNcols(R)) {
      bad_error("bad subscripts or zero matrix in Imatrix_submat()");
    } else {
      ImatrixMrows(R) = r;
      ImatrixMcols(R) = c;
    }
    return R;
}


/*
   **  Imatrix_print(M):  print an Imatrix
   **    if M is null print <<>> and return fail.
   **    otherwise print matrix and return true.
 */
Imatrix Imatrix_fprint(FILE *fout,Imatrix M)
{
    int i, j;

    if (M == 0) {
 #ifdef LOG_PRINT
       fprintf(fout,"<>\n")
#endif
;
        return 0;
    }
 #ifdef LOG_PRINT
   fprintf(fout,"<")
#endif
;
    for (i = 1; i <= ImatrixMrows(M); i++) {
        for (j = 1; j <= ImatrixMcols(M); j++) {
   #ifdef LOG_PRINT
         fprintf(fout,"%d", ImatrixMref(M, i, j))
#endif
;
            if (j < ImatrixMcols(M))  
#ifdef LOG_PRINT
                fprintf(fout,", ")
#endif
;
        }
        if (i < ImatrixMrows(M))
#ifdef LOG_PRINT
            fprintf(fout,";\n")
#endif
;
    }
#ifdef LOG_PRINT
    fprintf(fout,">\n")
#endif
;
    return M;
}

/*
** matrix_add(M1,M2,&M3) -- Add two Imatrixes:
** if M1, and M2 are incompatable (or null) complain and return false.
**  if *M3 has too little storage (or is null) free *M3 if nescesary
**                                            and create new storage.
*/
Imatrix Imatrix_add(int i1,Imatrix M1,int i2, Imatrix M2, Imatrix M3)
{
    int i, j;
    Imatrix R = M3;

    if (M1 == 0 || M2 == 0 || ImatrixMrows(M1) != ImatrixMrows(M2) || ImatrixMcols(M1) != ImatrixMcols(M2)) {
#ifdef LOG_PRINT
      fprintf(stderr, "matrix_add: dimensions dont match\n")
#endif
;
      return 0;
    }
    R=Imatrix_resize(R, ImatrixMrows(M1), ImatrixMcols(M1));
    for (i = 1; i <= ImatrixMrows(M1); i++)
      for (j = 1; j <= ImatrixMcols(M1); j++)
          ImatrixMref(R, i, j) =i1* ImatrixMref(M1, i, j) + i2*ImatrixMref(M2, i, j);
    return R;
}


/*
   **Imatrix_mull(M1,M2,&M3) -- multiply two Imatrixes:
   ** if M1, and M2 are incompatable (or null) complain and return false.
   ** if *M3 has too little storage (or is null) free *M3 if nescesary
   **                                           and create new storage.
   ** NOT USED YET
 */

Imatrix Imatrix_mul(Imatrix M1, Imatrix M2, Imatrix M3)
{
    int i, j, k;
    Imatrix R;

    if (M1 == 0 || M2 == 0 || ImatrixMcols(M1) != ImatrixMrows(M2)) {
#ifdef LOG_PRINT  
fprintf(stderr, "Imatrix_mull: incompatible matrices\n")
#endif
;
      return 0;
    }
    if (M3!=M1&&M3!=M2) R=Imatrix_resize(M3, ImatrixMrows(M1), ImatrixMcols(M2));
      else R=Imatrix_new(ImatrixMrows(M1),ImatrixMcols(M2));
    for (i = 1; i <= ImatrixMrows(M1); i++)
      for (j = 1; j <= ImatrixMcols(M2); j++) {
          ImatrixMref(R, i, j) = 0;
          for (k = 1; k <= ImatrixMcols(M1); k++)
            ImatrixMref(R, i, j) += ImatrixMref(M1, i, k) * ImatrixMref(M2, k, j);
      }
    if (M3==M1&&M3==M2) Imatrix_free(M3);
    return R;
}

/*
   ** Imatrix_dot(M1,M2,&M3) -- calculate M1*Transpose(M2):                
   ** if M1, and M2 are incompatable (or null) complain and return false.
   ** if *M3 has too little storage (or is null) free *M3 if nescesary 
   **                                           and create new storage.   
 */
Imatrix Imatrix_dot(Imatrix M1, Imatrix M2, Imatrix M3)
{
    int i, j, k;

    if (M1 == 0 || M2 == 0 || ImatrixMcols(M1) != ImatrixMcols(M2)) {
#ifdef LOG_PRINT
      fprintf(stderr, "Imatrix_dot: incompatible matrices\n")
#endif
;
      return 0;
    }
    M3 = Imatrix_resize(M3, ImatrixMrows(M1), ImatrixMrows(M2));

    for (i = 1; i <= ImatrixMrows(M1); i++)
      for (j = 1; j <= ImatrixMrows(M2); j++) {
          ImatrixMref(M3, i, j) = 0;
          for (k = 1; k <= ImatrixMcols(M1); k++)
            ImatrixMref(M3, i, j) += ImatrixMref(M1, i, k) * ImatrixMref(M2, j, k);
      }
    return M3;
}

/*
   ** dot_Ivector(M1,M2) -- vector dot product (integer valued)
   **  calculate dot product of first rows of M1, and M2.
   ** if M1, and M2 have are incompatable(or null) complain and return 0.
   ** NOT USED YET
 */
int dot_Ivector(Imatrix M1, Imatrix M2)
{
    int k, d = 0;
    if (M1 == 0 || M2 == 0 || ImatrixMcols(M1) != ImatrixMcols(M2)) {
#ifdef LOG_PRINT
      fprintf(stderr, "dot_Ivector: incompatible or null vectors\n")
#endif
;
      return 0;
    }
    for (k = 1; k <= ImatrixMcols(M1); k++)
      d += ImatrixVref(M1, k) * ImatrixVref(M2, k);
    return d;
}

/*
   ** Imatrix_equal(M1,M2)  
   **     return true if matrices represented by M1 and M2 are equal.
   **     false otherwise.  
 */
int Imatrix_equal(Imatrix M1, Imatrix M2)
{
    int i, j;
    if (M1 == 0 && M2 == 0)
      return TRUE;
    if (M1 == 0 || M2 == 0 || (ImatrixMrows(M1) != ImatrixMrows(M2)) || (ImatrixMcols(M1) != ImatrixMcols(M2))) {
      printf("in Imatrix_equal comparison failing\n");
      return FALSE;
    }
    for (i = 1; i <= ImatrixMrows(M1); i++)
      for (j = 1; j <= ImatrixMcols(M1); j++)
          if (ImatrixMref(M1, i, j) != ImatrixMref(M2, i, j))
            return FALSE;

    return TRUE;
}


int Imatrix_order(Imatrix M1, Imatrix M2)
{
    int i, j, c = 0;
    if (M1 == 0 && M2 == 0)
      return TRUE;
    if (M1 == 0 || M2 == 0 || (ImatrixMrows(M1) != ImatrixMrows(M2)) || (ImatrixMcols(M1) != ImatrixMcols(M2))) {
      printf("in Imatrix_equal comparison failing\n");
      return FALSE;
    }
    for (i = 1; i <= ImatrixMrows(M1); i++)
      for (j = 1; j <= ImatrixMcols(M1); j++)
          if ((c = ImatrixMref(M1, i, j) - ImatrixMref(M2, i, j)) != 0)
            return c;
    return 0;
}



int Imatrix_gcd_reduce(Imatrix M)
{
    int i, j, g = 0, gcd(int, int);

    if (M == 0)
      return 1;
    for (i = 1; i <= ImatrixMrows(M); i++)
      for (j = 1; j <= ImatrixMcols(M); j++) {
          if (g == 0 && ImatrixMref(M, i, j) != 0)
            g = abs(ImatrixMref(M, i, j));
          else if (ImatrixMref(M, i, j) != 0)
            g = abs(gcd(g, ImatrixMref(M, i, j)));
      }

    for (i = 1; i <= ImatrixMrows(M); i++)
      for (j = 1; j <= ImatrixMcols(M); j++)
          ImatrixMref(M, i, j) = ImatrixMref(M, i, j) / g;

    return g;
}


int gcd(int a, int b)
{
    int c;
    if (b == 0)
      return a;
    c = a % b;
    if (c == 0)
      return b;
    else
      return gcd(b, c);
}

/*
** int Imatrix_rref(Imatrix N, int *det)
**   replaces N by its reduced row echelon form (upper triangular).
**   returns the rank of the  matrix, and puts the determinant
**   of the first full rank maximal minor in det.
**   ALGORITHM 9.1 of Algorithms for Computer Algebra,
**      Geddes,Czapor,Labahn
**  NOTE: I have actually changed this to reduce by the gcd when posible
**        to protect against over flow. (for the C-integer version).
*/
#define M(i,j) (ImatrixMref(N,i,j))

int Imatrix_rref(Imatrix N,int *det){
  int divisor=1,sign=1,r=1,p,k,j,i,g;

  for(k=1;k<=ImatrixMcols(N)&&r<=ImatrixMrows(N);k++){
    for(p=r; p<=ImatrixMrows(N)&&M(p,k)==0;p++); 
    if (p<=ImatrixMrows(N)) {
        for(j=k;j<=ImatrixMcols(N);j++){i=M(p,j); M(p,j)=M(r,j); M(r,j)=i;}
        if (r!=p) sign*=-1;
        for(i=r+1;i<=ImatrixMrows(N);i++){
          g=gcd(M(i,k),M(r,k)); 
          for(j=k+1;j<=ImatrixMcols(N);j++){
            M(i,j)=(M(r,k)*M(i,j)-M(r,j)*M(i,k))/divisor; 
          }
          M(i,k)=0;
        }
        divisor=M(r,k);
        r++;
    }
  }
  *det=sign*divisor;
  return r-1;
}
/*
** int Imatrix_backsolve(Imatrix N, Imatrix S)
**  finds a solution to a triangular system of equations with more
**  equtions then unknowns
*/
int Imatrix_backsolve(Imatrix N, Imatrix S){
 int p,j,g,k;


 if (N==0 || S==0 || IMcols(N)!=IMcols(S)||IMcols(N)<=IMrows(N)) 
     bad_error("bad dimensions in backsolve");
 
 /* 
 ** find first missing pivot (in row p;  p=IMrows(N)+1 if no missing pivots
 */
 for(p=1;p<=ImatrixMrows(N) && ImatrixMref(N,p,p)!=0;p++) /* null body*/;
 /* 
 ** now set up for a solution to the AX=0 problem with A a p-1xp matrix whoose
 ** first pxp minor is nonsingular (use simple minded back substitution)
 */
 for(j=1;j<=ImatrixMcols(N);j++) ImatrixVref(S,j)=0;
 ImatrixVref(S,p)=1;
 for (j=p-1;j>=1;j--){
   for(k=j+1;k<=p;k++) ImatrixVref(S,j)-=ImatrixMref(N,j,k)*ImatrixVref(S,k);
   g=gcd(ImatrixVref(S,j),ImatrixMref(N,j,j));
   ImatrixVref(S,j)=ImatrixVref(S,j)/g;
   g=(ImatrixMref(N,j,j)/g);
   if (ImatrixVref(S,j)!=0) for(k=j+1;k<=p;k++) ImatrixVref(S,k)*=g;
 }
return 1;
}
#undef M



/*
** Imatrix_hermite:
**    Input: S an n by m Imatrix.
**   Output: True 
**
**   Sidefects: U points to a unitary nxn matrix
**              S gets replaced by its hermite normal form (U*S)
**            
**  Method:
**    U = Inxn
**    c = 1
**    while (c<=n) do
**      find mc>=c such that abs(S(mc,c)) has minimum non-zero value
**      if (mc!=c) then  switch rows mc and c (in S, and U)
**      if (S(c,c)<0) then multiply row c by -1
**      if (S(c,c)!=0) then 
**        for i in c+1:n  add S(i,c)/S(c,c) times row c to row i;
**      end(if)
**      if S(i,c)==0 for all i in c+1 ... n set c=c+1
**    end(while)
**
*/

int Imatrix_hermite(Imatrix S, Imatrix U){
  int c=1,mc,i,j,m,n,done,sign;
  int t=0,mv=0;
  m=ImatrixMcols(S);
  n=ImatrixMrows(S);

  if (S==0 || U==0 || ImatrixMrows(U)!=n || ImatrixMrows(U)!=n)
     bad_error("Incompatable matrices in Imatrix_hermite");

  /* Initialize U to nxn identity */
  U=Imatrix_resize(U,n,n);
  for(i=1;i<=n;i++){
    for(j=1;j<=n;j++){
      if (i==j) ImatrixMref(U,i,j)=1;
      else ImatrixMref(U,i,j)=0;
    }     
  }

  while(c<=n){
    /* find minimum entry in col c */
    mv=abs(ImatrixMref(S,c,c));
    mc=c;
    for(i=c+1;i<=n;i++){
      t=abs(ImatrixMref(S,i,c)); 
      if(mv==0 || (mv>t && t!=0)){
           mv=t;
           mc=i;
      }
    }

    /* if nescesary pivot to put min in row c and multiply by+-1
       to ensure diagonal entry is positive */
    if (mc!=c||ImatrixMref(S,mc,c)<0){
      if (ImatrixMref(S,mc,c)<0) sign=-1;
      else sign=+1;
      for(j=c;j<=m;j++){
        t=ImatrixMref(S,mc,j);
        ImatrixMref(S,mc,j)=ImatrixMref(S,c,j);
        ImatrixMref(S,c,j)=sign*t;
      }
      for(j=1;j<=n;j++){
        t=ImatrixMref(U,mc,j);
        ImatrixMref(U,mc,j)=ImatrixMref(U,c,j);
        ImatrixMref(U,c,j)=sign*t;
      }
    }

    /* if collumb is not zero do a reduction step */
    done=TRUE;
    if (ImatrixMref(S,c,c)!=0){
      for(i=c+1;i<=n;i++){
        t=ImatrixMref(S,i,c)/ImatrixMref(S,c,c);
        for(j=c;j<=m;j++){
           ImatrixMref(S,i,j)-=t*ImatrixMref(S,c,j);
        }
        for(j=1;j<=n;j++){
           ImatrixMref(U,i,j)-=t*ImatrixMref(U,c,j);
        }
        if (ImatrixMref(S,i,c)!=0) done=FALSE;
      }
    }
    /* if all entrees of col c bellow row c are zero go tonext col */
    if (done==TRUE) c++;
  }
  return TRUE;
 }

Imatrix Imatrix_dup(Imatrix M,Imatrix N){
  Imatrix R;
  int i,j;
  if (M==0) return 0;
  if (M==N) bad_error("Imatrix_dup: source and destination equal");
  R=Imatrix_resize(N,ImatrixMrows(M),ImatrixMcols(M));
  for(i=1;i<=ImatrixMrows(M);i++)
     for(j=1;j<=ImatrixMcols(M);j++) ImatrixMref(R,i,j)=ImatrixMref(M,i,j);
  return R;
}


int Imatrix_is_zero(Imatrix M){
  int i,j;
  if (M==0) return TRUE;
  for(i=1;i<=IMrows(M);i++){
       for(j=1;j<=IMcols(M);j++){
         if (*IMref(M,i,j)!=0) return FALSE;
       }
  }
  return TRUE;
}

/* end Imatrix.c */

/**************************************************************************/
/********************** implementations from Lists.c **********************/
/**************************************************************************/

void bad_error(char *);
void mem_free(void *);

node list_push(node item,node *stack){
   *stack=Cons(item,*stack);
   return item;
} 
  
int list_length(node L){
  int ct=0;
  while(L!=0){
    ct++;
    L=Cdr(L);
  }
  return ct;
}

node list_pop(node *stack){
   node temp;
   temp=Car(*stack);
   *stack=Cdr(*stack);
   return temp;
}
  
node list_first(node list){
  return Car(list);
}

node list_rest(node list){
  return Cdr(list);
}

int list_empty(node list){
  if (list==0) return TRUE;
  else return FALSE;
}

/* 
** Inserts a node g onto a list L in order (according to a comparison
** function comp) -- IF There is allready and equivalent object on
** the list and uniq is TRUE it is not inserted, and a value of zero 
** is returned, otherwise it is inserted and a value of 1 is returned.
*/

int list_insert(node g, node * L, int (*comp) (node, node), int uniq)
{
  int flg = 1;
  node ptr = *L;
  LOCS(2);
  PUSH_LOC(*L);
  PUSH_LOC(g);

  if ((*L == 0) || ((flg = comp(g, Car(*L))) >= 0)) {
    if (uniq==FALSE || flg != 0) *L = Cons(g, *L);
  } 
  else {
    while ((Cdr(ptr) != 0) &&
      ((flg = comp(g, Car(Cdr(ptr))))< 0)) {
       ptr = (node) Cdr(ptr);
    }
    if (uniq==FALSE || flg != 0)
      node_set_ptr(ptr,Cons(g, Cdr(ptr)), NODE, RIGHT);
  }
  POP_LOCS();
  return flg;
}

/* 
** list remove -- remove a specific elt from a list 
*/

int list_Imatrix_comp(node g1, node g2)
{
    if (g1 == 0 || g2 == 0 ||
    node_get_type(g1, LEFT) != IMTX || node_get_type(g2, LEFT) != IMTX) {
      bad_error("Non-IMTX nodes in sexpr_IMTX_comp()");
    }
    return Imatrix_order((struct Imatrix_t *) Car(g1),
                   (struct Imatrix_t *) Car(g2));
}


node list_cat(node l1,node l2){
  node ptr=l1;
  if (ptr==0) return l2;
  while(Cdr(ptr)!=0) ptr=Cdr(ptr);
  node_set_ptr(ptr,(void *)l2,NODE,RIGHT);
  return l1;
}



void xpl_fprint(FILE *fout,node L){
 node ptr;
 Dmatrix P;
 int i;
#ifdef LOG_PRINT
 fprintf(fout, "%% Solution list :\n");
 fprintf(fout, 
 "%%      Re(H) ,       Im(H) ,      Re(Xi) ,      Im(Xi) ,          T  \n");
 fprintf(fout, "%%\n");
 fprintf(fout,"{\n")
 #endif
;
ptr=L;
 while(ptr!=0){
   P=(Dmatrix)(Car(Car(ptr)));
 #ifdef LOG_PRINT
  fprintf(fout,"<");
   fprintf(fout,"%12g, %12g,",DVref(P,1),DVref(P,2))
#endif
;
   for(i=3;i<DVlength(P);i++){
#ifdef LOG_PRINT
         fprintf(fout," %12g,",DVref(P,i))
#endif
; 
         if (i%2==0 && i<DVlength(P)-1)
      #ifdef LOG_PRINT
           fprintf(fout,"\n                            ")

#endif
;
   }
 #ifdef LOG_PRINT
  fprintf(fout," %12g >",DVref(P,i))

#endif

;
   ptr=Cdr(ptr);
   if (ptr!=0) 
#ifdef LOG_PRINT
fprintf(fout,",")

#endif
;
 #ifdef LOG_PRINT
  fprintf(fout,"\n")
#endif
;
 }
#ifdef LOG_PRINT
 fprintf(fout,"};")
#endif
;
}

/* end Lists.c */

/**************************************************************************/
/********************** implementations from Pconfig.c ********************/
/**************************************************************************/

/*
** Pconfig.c                                Birk Huber 2-1-1994
**
**  Module for points and point configurations in affine space.
**
**  A point are represented by a node with a label(cstring) on the
**     left and a vector of coordinates(Imatrix) on the right.
**    
**       -----------------------
**       | PNT      | IMTX     |
**       -----------------------
**       |(char *)  | Imatrix  | ------> structure holding coords
**       -----------------------
**           |------------------------> string
** 
** A point configuration is represented by a node with a list of 
**      points on the right and the number of points in the list 
**      stored on the 
**      left.
**
**       -----------------------
**       | PCFG     |  NODE    |    
**       -----------------------        ------------------
**       | int      | Imatrix  | ------>| NODE   | NODE  | 
**       -----------------------        ------------------
**                                      | node   | node  | ----->
**                                      ------------------
**                                         |
**                                        point 1
 */

node pcfg_start(node P){ return Cdr(P);}
node pcfg_next_pnt(node *ptc){ 
    node pt;
    pt=Car(*ptc);
    *ptc=Cdr(*ptc);
    return pt;}

/*
** node pnt_new(char *s,Imatrix m)
**       input:  a string s 
**               an integer matrix m
**      output: a point with lable s, and coordinates m
*/
node pnt_new(char *s, Imatrix m)
{
    node R;
    R = node_new();
    node_set_ptr(R, s, PNT, LEFT);
    node_set_ptr(R, m, IMTX, RIGHT);
    return R;
}
/*
** pnt_free
**    input: a node
**   output: none
** effects: frees space allocated for lable and coordinates.
**          resets fields of node to null nodeptrs.
** error :  checks that n is a non-null pnt.
*/
void pnt_free(node n)
{
    if (n==0||node_get_type(n, LEFT) == PNT) {
        mem_free((char *) node_get_ptr(n, LEFT));
        Imatrix_free((Imatrix) node_get_ptr(n, RIGHT));
        node_set_ptr(n, (void *) 0, NODE, RIGHT);
        node_set_ptr(n, (void *) 0, NODE, LEFT);
    } else
        bad_error("error: pnt_free() called on non-point\n");
}
/*
** default display for points.
*/
int pnt_print(node n)
{
    if (n==0||node_get_type(n, LEFT) == PNT) {
      printf("(");
      printf("%s", (char *) node_get_ptr(n, LEFT));
      printf(",");
      Imatrix_print((Imatrix) node_get_ptr(n, RIGHT));
      printf(")\n");
      return 0;
    }
    bad_error("error: pnt_print() called on non-point\n");
    return 1;                 /*bad error causes abort */
}

/*
** pnt_is_point  
**       input: a node n
**       output:  TRUE if n is a point (and non-null)
**               FALSE  otherwise
*/
int pnt_is_point(node n){ 
      if (n!=0&&node_get_type(n,LEFT)==PNT) return TRUE;
      else return FALSE;
}
/*
** pnt_comp_lbl
**       input: two points;
**      output: value of strcmp on the lables of pt1 and pt2.
**      error conditions: checks that g1, and g2 are (non-null) points
*/
int pnt_comp_lbl(node g1, node g2){
  if ((pnt_is_point(g1)!=TRUE)||(pnt_is_point(g2)!=TRUE)){
    bad_error("Non-PNT nodes in pnt_comp_lbl()");
  }
  return strcmp(pnt_lable(g2),pnt_lable(g1));
}

/* 
** pnt_comp_rand
**    input:  two points
**   output:  1 or -1 randomly chosen.
**   error conditions: checks that g1 and g2 are (non-null)points
**
** NOTE: this is used to randomize point lists
*/
int pnt_comp_rand(node g1, node g2){
    if ((pnt_is_point(g1)!=TRUE)||(pnt_is_point(g2)!=TRUE)){
    bad_error("Non-PNT nodes in pnt_comp_rand()");
  }
  if (rand_int(0,1)==0) return -1;
  else return 1;
}

int pnt_comp_const(node g1, node g2){
  if ((pnt_is_point(g1)!=TRUE)||(pnt_is_point(g2)!=TRUE)){
    bad_error("Non-PNT nodes in pnt_comp_rand()");
  }
  return 1;
}

/*
** pnt_lable(node n);
**         input: a point n
**        outpur: the lable associated with n
**         error: checks that n is a (non-null) point
*/
char *pnt_lable(node n)
{
    if (n == 0 || node_get_type(n, LEFT) != PNT) {
      bad_error("error: pnt_label() called on non-point\n");
      /* will cause an exit */
    }
    return (char *) node_get_ptr(n, LEFT);
}

/*
** POINT CONFIGURATIONS
*/
/* 
** pcfg_new
**      input: none
**     output: an empty point configuration.
*/
node pcfg_new()
{
    node g;
    g = node_new();
    node_set_int(g, 0, PCFG, LEFT);
    node_set_ptr(g, (void *) 0, NODE, RIGHT);
    return g;
}

node pcfg_print(node n)
{
    if (n != 0 && node_get_type(n, LEFT) == PCFG) {
      printf("<");
      node_print((node)node_get_ptr(n, RIGHT));
      printf(">");
    } else
      bad_error("error: pcfg_print() called on non-PCFG\n");
    return n;
}

node pcfg_print_short(node n)
{
    if (n == 0 || node_get_type(n, LEFT) != PCFG)
      bad_error("error: pcfg_print_short() called on non-PCFG\n");
    n = (node) node_get_ptr(n, RIGHT);
    printf("<");
    while (n != 0) {
      printf("%s ", (char *) node_get_ptr(Car(n), LEFT));
      n = (node) node_get_ptr(n, RIGHT);
    }
    printf(">");
    return n;
}

/*
** pcfg_add
**   input:  a point "point"
**           a point config "config"
**   output: TRUE on successfull completion
**           (right now failure causes the program to abort)
**  side effects: point is added into config 
**                (which is kept in order according to point lables)
*/
int pcfg_add(node point, node config)
{   node ptr;
    LOCS(1);
    PUSH_LOC(config);
    if (config == 0 || 
         point == 0 || 
         node_get_type(config, LEFT) != PCFG ||
       node_get_type(point, LEFT) != PNT){
          bad_error("error: pcfg_add() called on non-PCFG\n");
    }
    ptr=Cdr(config);
    list_insert(point,&ptr,&(pnt_comp_lbl),FALSE);
    node_set_ptr(config,ptr,NODE,RIGHT);
    node_set_int(config, node_get_int(config, LEFT) + 1, PCFG, LEFT);
    POP_LOCS();
    return TRUE;
}
/*
** pcfg_remove
**   input:  a point "point"
**           a point configuration "config"
**           a pointer ptr locating "point" in "config" i.e.
**   output: TRUE if point is in config  
**           FALSE if point is not in config
**  side effect: if point is in config it is removed from config.
**  error conditions: checks that config is non-null pcfg
**                                pnt is non-null pnt
**  Note: ptr should point to the cons node for the point before
**        "point" in config's point list, if it doesn't then 
**        pcfg_remove will find the previous point itself
**           
*/
int pcfg_remove(node point, node config, node ptr)
{
  if (config == 0 ||point == 0||
        node_get_type(config, LEFT) != PCFG ||
        node_get_type(point, LEFT) != PNT)
     bad_error("error: pcfg_add() called on non-PCFG\n");
 /* check if location given for point is really right*/
  if (Car(Cdr(ptr))!=point) ptr=0;  
 /* find location of point if it is not allready known*/
  if (ptr==0)
     for(ptr=Cdr(config); Car(ptr)!=point && ptr!=0; ptr=Cdr(ptr)){
        ; /* null body*/
     }
  if (ptr==0) return FALSE; 
  node_set_ptr(ptr,Cdr(Cdr(ptr)),NODE,RIGHT);
  node_set_int(config,node_get_int(config,LEFT)-1,PCFG,LEFT);
  return TRUE;
}

/*
**pcfg_in(node point, node config);
**  test if point is in config. point must be given by a pointer
**  to a point in config. 
**
** Warning: Point must be given by a pointer to a point in config.
**          Point_config_in DOES NOT TELL wheather their is a point 
**          in Config with same coords/lable.
*/
int pcfg_in(node point, node config)
{
    if (config == 0 || point == 0 || node_get_type(config, LEFT) != PCFG ||
      node_get_type(point, LEFT) != PNT)
      bad_error("error: pcfg_in() called on non-PCFG\n");
    while ((config = Cdr(config)) != 0)
      if (point == Car(config))
          return TRUE;
    return FALSE;
}

/* 
** Imatrix pcfg_coords(node n, Imatrix R)
**    returns a matrix whoose rows are the coordinates of
**    points in the configuration. Space for the matrix is
**    allocated from R if possible, otherwise R is freed and
**    space realocated.
*/

Imatrix pcfg_coords(node n, Imatrix R)
{
    int i, j;
    node ptr = n;
    if (n == 0 || node_get_type(n, LEFT) != PCFG)
      bad_error("error: pcfg_coords() called on non-PCFG\n");
    R = Imatrix_resize(R, pcfg_npts(n), pcfg_dim(n));
    for (i = 1; i <= pcfg_npts(n); i++) {
      ptr = Cdr(ptr);
      for (j = 1; j <= pcfg_dim(n); j++)
          *(IMref(R, i, j)) = *(IVref(pnt_coords(Car(ptr)), j));
    }
    return R;
}

/* 
** Imatrix pcfg_M(node n, Imatrix R)
**    returns a matrix whoose rows are the coordinates of
**    points in the configuration, after translation to put the 
**    first point at the origen.    Space for the matrix is
**    allocated from R if possible, otherwise R is freed and
**    space realocated.
*/
Imatrix pcfg_M(node n, Imatrix R)
{
    int i, j;
    Imatrix P0;
    node ptr = n;



    if (n == 0 || node_get_type(n, LEFT) != PCFG)
      bad_error("error: pcfg_M called on non-PCFG\n");
    R = Imatrix_resize(R, pcfg_npts(n) - 1, pcfg_dim(n));
    ptr = Cdr(ptr);
    P0 = pnt_coords(Car(ptr));
    for (i = 1; i <= pcfg_npts(n) - 1; i++) {
      ptr = Cdr(ptr);
      for (j = 1; j <= pcfg_dim(n); j++)
          *(IMref(R, i, j)) = *(IVref(pnt_coords(Car(ptr)), j))
            - *(IVref(P0, j));
    }
    return R;
}
/* 
** node pcfg_face(node PC, Imatrix norm);
**   returns the point configuration consisting of points 
**   of PC which lie on the face suported by N.
*/
node pcfg_face(node PC, Imatrix norm)
{
    node Face = 0, ptr = 0;
    Imatrix M = 0, D = 0;
    int s, i;
    LOCS(2);

    /* save locals */
    PUSH_LOC(PC);
    PUSH_LOC(Face);
    M = pcfg_coords(PC, M);
    D = Imatrix_dot(norm, M, D);
    s = *IMref(D, 1, 1);
    for (i = 2; i <= IMcols(D); i++)
      if (*IMref(D, 1, i) < s)
          s = (*IMref(D, 1, i));
    Face = pcfg_new();
    ptr = PC;
    for (i = 1; i <= IMcols(D); i++) {
      ptr = Cdr(ptr);
      if (*IMref(D, 1, i) == s)
          pcfg_add(Car(ptr), Face);
    }
    Imatrix_free(D);
    Imatrix_free(M);
    POP_LOCS();
    return Face;
}
/*
**  is_normal_good(Imatrix normal, Imatrix N)
**  determine if the (inner normal) associated to a facet
**  satasfies the sign conditions specified in N.
**
**  each coordinate of the (inner)normal satasfy the condition
**  specified by the corresponding coordinate of N:
**           1:   normal[i]  must be non-negative
**           2:   normal[i]  must be positive
**          -1:   normal[i]  must be non-positive
**          -2:   normal[i]  must be negative
**           0:   no condition required of normal[i]
*/

int is_normal_good(Imatrix normal, Imatrix N)
{
    int i;
    if (N == 0)
      return TRUE;
    for (i = 1; i <= IVlength(N); i++) {
      if (*(IVref(N, i)) == 0)      /* Skip */
          ;
      else if (*(IVref(N, i)) * (*IVref(normal, i)) < 0)
          return FALSE;
      else if ((*(IVref(normal, i)) == 0) && (abs(*IVref(N, i)) == 2))
          return FALSE;
    }
    return TRUE;
}

/* end Pconfig.c */

/**************************************************************************/
/******************** implementations from Pcomplex.c *********************/
/**************************************************************************/

#ifndef PI
#define PI (double)3.14159265358979323846264338328
#endif


fcomplex Cadd(fcomplex a,fcomplex b)
{     fcomplex c;
      c.r=a.r+b.r;
      c.i=a.i+b.i;
      return c;
}

fcomplex Csub(fcomplex a,fcomplex b)
{     fcomplex c;
      c.r=a.r-b.r;
      c.i=a.i-b.i;
      return c;
}

fcomplex Cmul(fcomplex a,fcomplex b)
{     fcomplex c;
      c.r=a.r*b.r-a.i*b.i;
      c.i=a.i*b.r+a.r*b.i;
      return c;
}

              
fcomplex ItoC(int i)
{ fcomplex c;
  c.r=i;
  c.i=0.0;
  return c;
}

fcomplex DtoC(double i)
{ fcomplex c;
  c.r=i;
  c.i=0.0;
  return c;
}

fcomplex Complex(double re, double im)
{     fcomplex c;
      c.r=re;
      c.i=im;
      return c;
}

fcomplex Conjg(fcomplex z)
{     fcomplex c;
      c.r=z.r;
      c.i = -z.i;
      return c;
}

fcomplex Cdiv(fcomplex a,fcomplex b)
{     fcomplex c;
      double r,den;
      if (fabs(b.r) >= fabs(b.i)) {
            r=b.i/b.r;
            den=b.r+r*b.i;
            c.r=(a.r+r*a.i)/den;
            c.i=(a.i-r*a.r)/den;
      } else {
            r=b.r/b.i;
            den=b.i+r*b.r;
            c.r=(a.r*r+a.i)/den;
            c.i=(a.i*r-a.r)/den;
      }
      return c;
}

fcomplex Cpow(fcomplex a, int d)
  { int i;
    fcomplex c;
    c.r=1.0;
    c.i=0.0;
    if (d<0) { a=Cdiv(c,a);
               d*=-1;}
    for (i=1;i<=d;i++) c=Cmul(c,a);
    return c;
    }
double Cabs(fcomplex z)
{     double x,y,ans,temp;
      x=fabs(z.r);
      y=fabs(z.i);
      if (x == 0.0)
            ans=y;
      else if (y == 0.0)
            ans=x;
      else if (x > y) {
            temp=y/x;
            ans=x*sqrt(1.0+temp*temp);
      } else {
            temp=x/y;
            ans=y*sqrt(1.0+temp*temp);
      }
      return ans;
}

fcomplex Csqrt(fcomplex z)
{     fcomplex c;
      double x,y,w,r;
      if ((z.r == 0.0) && (z.i == 0.0)) {
            c.r=0.0;
            c.i=0.0;
            return c;
      } else {
            x=fabs(z.r);
            y=fabs(z.i);
            if (x >= y) {
                  r=y/x;
                  w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
            } else {
                  r=x/y;
                  w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
            }
            if (z.r >= 0.0) {
                  c.r=w;
                  c.i=z.i/(2.0*w);
            } else {
                  c.i=(z.i >= 0) ? w : -w;
                  c.r=z.i/(2.0*c.i);
            }
            return c;
      }
}

fcomplex RCmul(double x, fcomplex a)
{     fcomplex c;
      c.r=x*a.r;
      c.i=x*a.i;
      return c;
}

void printC(fcomplex c)
{     if (c.r!=0 && c.i !=0) {
            if (c.i>0) printf("(%g+%g*I)",c.r,c.i);
            else  printf("(%g-%g*I)",c.r,fabs(c.i));
      }
      else if (c.r!=0) printf("%g",c.r);
      else if (c.i!=0) printf("(%g*I)",c.i);
        else printf("0.0");
}

/*-------------------------------------------------------------------------
function  (fcomplex) Croot
  arguments:   n an integer;
               z a complex number;
  returns:     a primitive nth root of z;
-----------------------------------------------------------------------*/


fcomplex RootOfOne(int j, int n)
{ fcomplex C;
  C=Complex(cos((2.0*PI*j)/n),sin((2.0*PI*j)/n));
  return C;}



/*-------------------------------------------------------------------------
function  (fcomplex) Croot
  arguments:   n an integer;
               z a complex number;
  returns:     a primitive nth root of z;
-----------------------------------------------------------------------*/
fcomplex Croot(fcomplex z, int n)
  {   fcomplex w;
      double mod,arg,ni=1.0/n;

      if ((z.i==0)&&(z.r==0)) w=Complex(0.0,0.0);
         else{ mod=Cabs(z);
               if (fabs(z.i)>=fabs(z.r)){
                      arg=acos(z.r/mod);
                      if (z.i<0) arg=-1*arg;}
                else { arg=asin(z.i/mod);
                      if (z.r<0) arg=PI-arg;}
               arg=arg/n;
               mod=fabs(pow(mod,ni));
               w=Complex(cos(arg),sin(arg));
               w=RCmul(mod,w); }

     return w;}

/* end Pcomplex.c */

/**************************************************************************/
/*********************** implementations from Poly.c **********************/
/**************************************************************************/

void bad_error(char *);   

int *poly_exp(monomial m, int i){
  if (m==0 || i<1 || m->R->n < i ) {
      printf("i=%d, n=%d\n",i,m->R->n);
      bad_error("index out of bounds in monomial");
 }
  return &(m->exps[i-1]);
}   

int poly_deg(polynomial1 p){
    int i,n,deg=0,tdeg;
    n=p->R->n;
    for(i=0;i<n;i++) deg+=p->exps[i];
    while((p=p->next)!=0){
      tdeg=0;
      for(i=0;i<n;i++) tdeg+=p->exps[i];
      if (tdeg>deg) deg=tdeg;
    }
   return deg;
}

int *poly_homog(polynomial1 p){
 return &(p->homog);
}

void ring_set_var(Pring R, int n, char *lable){
   strncpy(R->vars[n],lable,RING_VAR_L);
}
char *ring_var(Pring R,int i){return R->vars[i];}
void ring_set_def(Pring R,  char *lable){
   strncpy(R->def,lable,RING_VAR_L);
}
char *ring_def(Pring R){return R->def;}
int poly_dim(monomial m){  
  if (m==0) bad_error("null monomial in poly_def");
  return m->R->n;
}
int ring_dim(Pring R){  
  if (R==0) bad_error("null Ring in ring_dim");
  return R->n;
}
int *poly_def(monomial m){  
  if (m==0) bad_error("null monomial in poly_def");
  return &(m->def);
}
fcomplex *poly_coef(monomial m){
  if (m==0) bad_error("null monomial in poly_coef");
  return &(m->coef);
}
 monomial poly_next(monomial m){
  if (m==0) bad_error("null monomial in poly_next");
  return m->next;
}
monomial poly_set_next(monomial m,monomial m2){
  if (m==0) bad_error("null monomial in poly_next");
  return (m->next=m2);
}

Pring poly_ring(monomial m){
  if (m==0) bad_error("null monomial in poly_next");
  return m->R;
}


Pring makePR(int n)

{ 
      int i;
      Pring a;
      a=(Pring)mem_malloc(sizeof(struct Pring_tag)); 
      if (a==0) bad_error("malloc failed 1 in make_PR");
      a->n=n;
      a->vars=(char **)mem_malloc(n*sizeof(char *));  
      if (a->vars==0) bad_error("malloc failed 2 in make_PR");
      for(i=0;i<n;i++) { 
            a->vars[i]=(char *)mem_malloc(RING_VAR_L*sizeof(char)); 
            if (a->vars==0) bad_error("malloc failed 3 in make_PR");
      }
        a->def=(char *)mem_malloc(RING_VAR_L*sizeof(char)); 
        if (a->vars==0) bad_error("malloc failed 4 in make_PR");
    return a;
}

Pring free_Pring(Pring R)
 {
 int i;
 mem_free(R->def);  
 for (i=0;i<R->n;i++){ mem_free(R->vars[i]);}
 mem_free( (char *) R->vars);
 mem_free( R); 
 return 0;
 }
       
polynomial1 makeP(Pring R)

{ 
  polynomial1 m;
  int i,n;

  n=R->n;
  m=(polynomial1)mem_malloc(sizeof(struct mono_tag));
  if (m==0){ printf("malloc 1 failed: makeP\n"); exit(2); }
  m->R=R;
  m->exps=(int *)mem_malloc(n*sizeof(int));
  if (m->exps==0){printf(" malloc 2 failed: makeP\n"); exit(2);}
  for(i=0;i<n;i++) m->exps[i]=0;
  m->coef=Complex(0.0,0.0);
  m->def=0;
  m->next=0;
  return m;
}
     
polynomial1 freeP(polynomial1 p)
{
  polynomial1 m=p;
  while (m!=0){
      p=p->next;
      mem_free((char *)m->exps); 
      mem_free((char *)m); 
      m=p;
  }
  return 0;
}

polynomial1 copyM(polynomial1 p)
{
  polynomial1 p1;
  int i;
  if (p==0) return 0;
  p1=makeP(p->R);
  p1->coef=p->coef;
  p1->def=p->def;
  for (i=0;i<p->R->n;i++) p1->exps[i]=p->exps[i];
 return p1;
}

 
polynomial1 copyP(polynomial1 p)
{   
  polynomial1 p1,p2;
  int i;
      
  if (p==0) return 0;
  p1=makeP(p->R);
  p2=p1;
  while (p!=0){
    if (p->next!=0) p2->next=makeP(p->R);
    p2->coef=p->coef;
    p2->def=p->def;
    for (i=0;i<p->R->n;i++) p2->exps[i]=p->exps[i];
    p=p->next;
    p2=p2->next;
  }
 return p1;
}
    
void printP(polynomial1 P)
{ 
  int i,zt=0;

  /* DEBUG */
  /* printf(" In "); */
 
  while (P!=0){
     if (P->coef.i!=0 || P->coef.r!=0){
        printC(P->coef);
        zt=1;

      /*    printf("("); */

        for (i=0;i<P->R->n;i++) { 
      
        /*       printf("\n(P->exps([%d])= %d\n", i,P->exps[i]);    */
        /*
        printf("%d", P->exps[i]);
        if (i == P->R->n - 1)
          printf(")");
        else
          printf(",");
        */

       if (P->exps[i]!=0) {
           zt=1;
         
         printf("*%s",P->R->vars[i]);
         if (P->exps[i]!=1) {
           /*              printf("^{without P->def-loop}");   */
              if (P->exps[i]>0) 
            printf("^%d",P->exps[i]);
              else 
            printf("(%d)",P->exps[i]);
           }
        }
     }
     if (P->def!=0){
 

      printf("*%s\n",P->R->def);
        if (P->def !=1) {
             printf("^{with  P->def-loop}");
             if (P->def >0) printf("%d ",P->def);
              else printf("(%d){second} ",P->def);
             }
     } 
     /*     if (P->next != 0) printf("+{new mono}\n");   */
     if (P->next != 0) printf(" + ");
     }
     P=P->next;
  }
  if (zt==0) printf("0.0nana");

  /* DEBUG */
  /* printf("  Out  "); */
} 

int orderMM(polynomial1 P1,polynomial1 P2)
{
  int i,d;
  if (P1==0 && P2==0) return 0;
   else if (P1==0) return -1;
    else if (P2==0) return 1;
  if (P1->R != P2->R) 
    bad_error("error in order: monomials must belong to same ring");
      
  d=P1->def-P2->def;
  if (d>0) return 1;
  else if (d<0) return -1;

  for (i=0;i<P1->R->n;i++) { 
    d=P1->exps[i]-P2->exps[i];
    if (d>0) return 1;
     else if (d<0) return -1;
  }
  return 0;
}

int orderPP(polynomial1 P1,polynomial1 P2)
{
      int i,d;
      while (P1!=0 && P2!=0){
                 d=P1->def-P2->def;
                 if (d>0) return 1;
                 else if (d<0) return -1;

            for (i=0;i<P1->R->n;i++) {
                  d=P1->exps[i]-P2->exps[i];
                         if (d>0) return 1;
                  else if (d<0) return -1;
            }
                P1=P1->next;
                P2=P2->next;
      }
        if (P1!=0) return 1;
         else if (P2 !=0) return -1;
          else return 0;
}
   
polynomial1 ItoP(int c,Pring R)
{ 
  polynomial1 p;
  p=makeP(R);
  p->coef.r=c;
  return p;
}

polynomial1 DtoP(double c, Pring R)
{ 
  polynomial1 p;
  p=makeP(R);
  p->coef.r=c;
  return p;
}

polynomial1 CtoP(fcomplex c, Pring R)
{ 
  polynomial1 p;
  p=makeP(R);
  p->coef=c;
  return p;
}
 
polynomial1 addPPP(polynomial1 P1, polynomial1 P2, polynomial1 P3)
{
  polynomial1 pt;
  struct mono_tag A;
  int d,i;
      
  pt=&(A);
if (P1==0){ if (P3==0) freeP(P3);
            return copyP(P2);
          }
if (P2==0){ if (P3==0) freeP(P3);
            return copyP(P1);
          }
/*
if (P1->R != P2->R) bad_error(" polynomial1s in addPPP must have equal rings");
*/ 
 A.R=P1->R;
  A.next=0;
  while(P1!=0&&P2!=0){
    d=orderMM(P1,P2);
    if (d==0) {
          A.coef=Cadd(P1->coef,P2->coef);
          if (Cabs(A.coef)!=0){
            pt->next=makeP(A.R);
            pt=pt->next;
            pt->coef=A.coef;
            pt->def=P1->def;
            for(i=0;i<A.R->n;i++)pt->exps[i]=P1->exps[i];
          }
          P1=P1->next;
          P2=P2->next;
    }
    else if (d==-1){
          pt->next=makeP(A.R);
          pt=pt->next;
          pt->coef=P2->coef;
          for(i=0;i<A.R->n;i++)pt->exps[i]=P2->exps[i];
          pt->def=P2->def;
          P2=P2->next;
          }
    else {
          pt->next=makeP(A.R);
        pt=pt->next;
          pt->def = P1->def;
          pt->coef=P1->coef;
          for(i=0;i<A.R->n;i++)pt->exps[i]=P1->exps[i];
          P1=P1->next;
      }
  }

if (P1!=0) pt->next=copyP(P1);
 else if (P2!=0) pt->next=copyP(P2);

if (P3!=0) freeP(P3);

if (A.next==0) A.next=makeP(A.R);
return A.next;
}           

polynomial1 subPPP(polynomial1 P1, polynomial1 P2, polynomial1 P3)
{
  polynomial1 pt;
  struct mono_tag A;
  int d,i;

  pt=&(A);
if (P1!=0 || P2!=0){
  if (P1->R != P2->R) bad_error(" polynomial1s in addPPP must have equal rings");
  A.R=P1->R;
  A.next=0;
  while(P1!=0&&P2!=0){
    d=orderMM(P1,P2);
    if (d==0) {
          A.coef=Csub(P1->coef,P2->coef);
          if (Cabs(A.coef)!=0){
            pt->next=makeP(A.R);
            pt=pt->next;
            pt->coef=A.coef;
            pt->def=P1->def;
            for(i=0;i<A.R->n;i++)pt->exps[i]=P1->exps[i];
          }
          P1=P1->next;
          P2=P2->next;
    }
    else if (d==-1){
          pt->next=makeP(A.R);
          pt=pt->next;
          pt->coef=RCmul(-1.0,P2->coef);
          pt->def=P2->def;
          for(i=0;i<A.R->n;i++)pt->exps[i]=P2->exps[i];
          P2=P2->next;
          }
    else {
          pt->next=makeP(A.R);
          pt=pt->next;
          pt->coef=RCmul(1.0,P1->coef);
          pt->def=P1->def;
          for(i=0;i<A.R->n;i++)pt->exps[i]=P1->exps[i];
          P1=P1->next;
        }
  }
}
if (P1!=0) pt->next=copyP(P1);
 else if (P2!=0) pt->next=mulCPP(Complex(-1.0,0.0),P2,0);

if (P3!=0) freeP(P3);

if (A.next==0) A.next=makeP(A.R);
return A.next;
}

polynomial1 mulCPP(fcomplex c, polynomial1 P1, polynomial1 P2)
{
  polynomial1 p;

  if (P1==0) { if (P2!=0) freeP(P2);
                          return 0;}
  if (P2!=P1 ) { if ( P2 !=0)  freeP(P2); 
                 P2=copyP(P1);}
  p=P2;
  while (p!=0){
    p->coef=Cmul(c,p->coef);
    p=p->next;
  }
 return P2;
}

polynomial1 divCPP(fcomplex c, polynomial1 P1, polynomial1 P2)
{
  polynomial1 p;

  if (P1==0) { if (P2!=0) freeP(P2);
                          return 0;}
  if (P2!=P1 ) { if ( P2 !=0)  freeP(P2);
                 P2=copyP(P1);}
  p=P2;
  while (p!=0){
    p->coef=Cdiv(p->coef,c);
    p=p->next;
  }
 return P2;
}


polynomial1 mulMPP(polynomial1 mi, polynomial1 P1, polynomial1 P2)
{
  polynomial1 p,m=mi;
  int i,free_m=0;

m->R=P1->R;

if (P1==0||m==0) { if (P2!=0) freeP(P2);
                      return 0;}
  if (P1->R!=m->R){ printf("\nRings must be equal in mulMPP");
                   if (P2!=0) freeP(P2);
                          return 0;}
  if (P2==m) free_m = 1;
  if (P2!=P1 ) { if ( P2 !=0 && P2 !=m)  freeP(P2); 
                 P2=copyP(P1);}
  p=P2;
  while (p!=0){
    p->coef=Cmul(m->coef,p->coef);
    p->def+=m->def;
    for(i=0;i<m->R->n;i++) p->exps[i]+=m->exps[i];
    p=p->next;
  }
 if (free_m==1) freeP(m);
 return P2;
}

polynomial1 divMPP(polynomial1 mi, polynomial1 P1, polynomial1 P2)
{
  polynomial1 p,m=mi;
  int i,free_m=0;

  if (P1==0||m==0) { if (P2!=0) freeP(P2);
                      return 0;}
  if (P1->R!=m->R){ printf("Rings must be equal in divMPP");
                   if (P2!=0) freeP(P2);
                          return 0;}
  if (m->next!=0){ printf("divisor must be a monomial in divMPP");
                   if (P2!=0) freeP(P2);
                    return 0;}

  if (P2==m) free_m = 1;
  if (P2!=P1 ) { if ( P2 !=0 && P2 !=m)  freeP(P2); 
                 P2=copyP(P1);}
  p=P2;
  while (p!=0){
    p->coef=Cdiv(p->coef,m->coef);
    p->def-=m->def;
    for(i=0;i<m->R->n;i++) p->exps[i]-=m->exps[i];
    p=p->next;
  }
 if (free_m==1) freeP(m);
 return P2;
}


polynomial1 mulPPP(polynomial1 P1, polynomial1 P2, polynomial1 P3)
{
  polynomial1 pt2,pt3=0;
    
  if (P1==0 || P2==0) { 
    if (P3!=0) freeP(P3);
    return 0;
  }
  /*
  if (P1->R != P2->R) {
    if (P3!=0) freeP(P3);
    bad_error("error in addPP: unequal rings");
  }
  */
  while(P1!=0) { 
    pt2=mulMPP(P1,P2,0);
    pt3=addPPP(pt2,pt3,pt3);
    freeP(pt2);               /* this realy should be done more efficiently */
    P1=P1->next;
    }

 if (P3!=0) free(P3);
 return pt3;
}


                  
polynomial1 expIPP(int x, polynomial1 P, polynomial1 P3)
{
  polynomial1 m,tmp;
      int i;
      if (P==0) bad_error("undefined poly in powPIP");
      if (x==0) {m=makeP(P->R);
                 m->coef.r=1.0;
                 if (P3!=0) freeP(P3);
                 return m;}
      if (P->next==0) { m=copyP(P);
                       m->coef=Cpow(m->coef,x);
                       m->def*=x;
                       for(i=0;i<m->R->n;i++) m->exps[i]*=x;
                       if(P3!=0) freeP(P3);
                       return m;}
      if (x<0) bad_error("powPIP can not raise a polynomial1 to a negative power");
      m=copyP(P);
      for(i=2;i<=x;i++) {
                       tmp=mulPPP(m,P,0);
                       freeP(m);
                       m=tmp;}
       if (P3!=0) freeP(P3);
       return m;
       }
       
polynomial1 unliftP(polynomial1 p)
{
  polynomial1 p1,p2;
  int i;

  if (p==0) return 0;
  p1=makeP(p->R);
  p2=p1;
  while (p!=0){
    if (p->next!=0) p2->next=makeP(p->R);
    p2->coef=p->coef;
    p2->def=0;
    for (i=0;i<p->R->n;i++) p2->exps[i]=p->exps[i];
    p=p->next;
    p2=p2->next;
  }
 return p1;
}

polynomial1 Homogenize(polynomial1 pi,Pring R)
{ 
  polynomial1 pt,p;
  int i = 0;
  int d = 0;
  int dp = 0;

  if (pi==0) bad_error("Homogenize called on null polynomial1");
  p=copyP(pi);
  pt=p;
  while (pt!=0){
    d=0;
    for (i=0;i<p->R->n-1;i++) d+=pt->exps[i];
    if (p==pt) dp=d;
      else if (dp<d) dp=d;
    pt->R=R;
    pt=pt->next;
  }
  pt=p;
  while (pt!=0){
    d=0;
    for (i=0;i<p->R->n-1;i++) d+=pt->exps[i];        
    pt->exps[p->R->n-1]=dp-d;
    pt=pt->next;
  }
return p;}

polynomial1 Prog_Eq(Pring R){
 polynomial1 p,pt;
 int i;
 pt=(p=ItoP(-1,R));
 for(i=0;i<R->n;i++) {
   pt->next=ItoP(1,R);
   pt=pt->next;
   pt->exps[i]=2;
 }
 return p;
}

/* end Poly.c */

/**************************************************************************/
/********************** implementations from utime.c **********************/
/**************************************************************************/

int set_mark(){
  return clock();
}

int read_mark(int timeset){
 return (clock()-timeset) / CLOCKS_PER_SEC;
}

/* end utime.c */

/**************************************************************************/
/*********************** implementations from Aset.c **********************/
/**************************************************************************/

/* 
** An Aset is an r-tupple of point configurations. It represents
** the support set of a polynomial system, and is a presentation for
** the minkowski sum of the convex hulls of its constituent
** point configurations. 
** 
** Asets are constructed from nodes as follows:
** _____________     ____________     _____________       ------------
**| ASET | NODE |   | INT | NODE |   | NODE | NODE |     | NODE | NODE
** -------------     ------------     -------------       ------------
**|  R   | node-|-->|  D  | node |-->| node | node |->...| node | NULL
** -------------     ------------     --|----------       --|---------
**                                    point config 1.   point config R
**                                                       
**
*/

#define R(A) (node_get_int(A,LEFT))
#define D(A) ((int)node_get_int((node)node_get_ptr(A,RIGHT),LEFT))
#define St(A) ((node)(node_get_ptr(A,RIGHT)))
#define C(P) ((node)(node_get_ptr(P,LEFT)))

node aset_new(int r, int d)
{
    node A = 0, ptr = 0;
    LOCS(2);
    PUSH_LOC(A);
    PUSH_LOC(ptr);
    A = node_new();
    node_set_int(A, r, ASET, LEFT);
    node_set_ptr(A, node_new(), NODE, RIGHT);
    node_set_int((node) Cdr(A), d, INT, LEFT);
    while ((r--) > 0)
      ptr = Cons(pcfg_new(), ptr);
    node_set_ptr((node) Cdr(A), ptr, NODE, RIGHT);
    POP_LOCS();
    return A;
}

int aset_r(node A)
{
    return R(A);
}
int aset_dim(node A)
{
    return D(A);
}
int aset_npts(node A)
{
    int i = 0;
    A = St(A);
    while ((A = Cdr(A)) != 0)
      i += pcfg_npts(Car(A));
    return i;
}
node aset_start_cfg(node ptr)
{
    return Cdr(ptr);
}
node aset_next_cfg(node * ptr)
{
    *ptr = Cdr(*ptr);
    return Car(*ptr);
}

node aset_start_pnt(node ptr)
{
    return Cdr(ptr);
}

node aset_next_pnt(node * ptr)
{
    *ptr = Cdr(*ptr);
    return Car(*ptr);
}

int aset_add(node A, int R, node point)
{
    node ptr = A;
    LOCS(1);
    PUSH_LOC(A);
    if (A == 0)
      bad_error("Null Aset in aset_add()");
    if (point == 0)
      bad_error("Null point in aset_add()");
    if (R >= 1 && R(A) >= R) {
      ptr = St(A);
      while ((R--) > 0)
          ptr = Cdr(ptr);
      pcfg_add(point, Car(ptr));
    } else
      bad_error("index out of bounds in aset_add()");
    POP_LOCS();
    return 1;
}

node aset_new_pt(int N, char *lable)
{
    Imatrix C;
    node point;
    C = Ivector_new(N);
    point = pnt_new(lable, C);
    return point;
}

int aset_pnt_set(node point, int i, int j)
{
    return (*IVref((Imatrix) Cdr(point), i) = j);
}
int aset_pnt_get(node point, int i)
{
    return *IVref((Imatrix) Cdr(point), i);
}

int aset_unlift(node A)
{
    int d;
    node ptp, ptc;

    d = aset_dim(A);
    A = St(A);
    while ((ptc = aset_next_cfg(&A)) != 0)
      while ((ptp = aset_next_pnt(&ptc)) != 0)
          aset_pnt_set(ptp, d, 0);

    return 1;
}

int aset_randlift(node A, int seed, int L, int U)
{
    int d;
    node ptp, ptc;

    rand_seed(seed);
    d = aset_dim(A);
    A = St(A);
    while ((ptc = aset_next_cfg(&A)) != 0)
      while ((ptp = aset_next_pnt(&ptc)) != 0)
          aset_pnt_set(ptp, d, rand_int(L, U));

    return 1;
}

node aset_print(node A)
{
 int R,D,N,i,first=1;
 Imatrix Coords;
 node ptr,ptc;

 if (A != 0 && node_get_type(A, LEFT) == ASET) {
   R = R(A);
   D = D(A);
   ptr = St(A);
   while ((ptr = (node) Cdr(ptr)) != 0) {
     first=1;
     ptc=C(ptr);
     N=pcfg_npts(ptc);
     while((ptc=Cdr(ptc))!=0){
        if (first==1){
             printf("   <");
             first=0;
        }
        else printf("    ");
        Coords=pnt_coords(Car(ptc));    
        printf("<");
        for(i=1;i<=D;i++) {
           printf(" %d",*IVref(Coords,i));
           if (i<D) printf(",");
        }
        if (--N>0) printf(" >, ");
        else {
          printf(" >>");
          if (--R > 0) printf(",");
          else printf(" ");
        }
        printf("       %% %s \n",(char *)pnt_lable(Car(ptc)));
     }
   }
 } else if (A != 0)
    bad_error("error: aset_print() called on non-ASET\n");
   return A;
}

node aset_print_short(node A)
{
 int R,N,first=1;  
 node ptr,ptc;
 if (A != 0 && node_get_type(A, LEFT) == ASET) {
   R = R(A);
   ptr = St(A);
   printf("{ ");
   while ((ptr = (node) Cdr(ptr)) != 0) {
     first=1;
     ptc=C(ptr);
     N=pcfg_npts(ptc);
     while((ptc=Cdr(ptc))!=0){
        if (first==1){
             printf("{");
             first=0;
        }
        printf(" %s",(char *)pnt_lable(Car(ptc)));   
        if (--N ==0) printf("}");
        else printf(",");
     }
     if (--R>0)printf(",");
     else printf(" }");
   }
 } else if (A != 0)
    bad_error("error: aset_print() called on non-ASET\n");
   return A;
}

node aset_face(node A, Imatrix N)
{
    node ptr = A, ptc = 0, res = 0, rptr = 0;
    LOCS(2);
    PUSH_LOC(A);
    PUSH_LOC(res);
    res = aset_new(aset_r(A), aset_dim(A));
    ptr = St(A);
    rptr = St(res);
    while ((ptc = aset_next_cfg(&ptr)) != 0) {
      rptr = Cdr(rptr);
      node_set_ptr(rptr, pcfg_face(ptc, N), NODE, LEFT);
    }
    POP_LOCS();
    return res;
}

Imatrix aset_type(node A, Imatrix T)
{
    int r = 1, i;
    Imatrix M = 0;
    node ptr, ptc;
    LOCS(1);
    PUSH_LOC(A);
    ptr = aset_start_cfg(A);
    T = Imatrix_resize(T, 1, aset_r(A));
    while ((ptc = aset_next_cfg(&ptr)) != 0) {
      M = pcfg_M(ptc, M);
      for (i = 1; i <= IMrows(M); i++)
          *IMref(M, i, IMcols(M)) = 0;
      *IVref(T, r++) = Imatrix_rref(M, &i);
    }
    Imatrix_free(M);
    POP_LOCS();
    return T;
}

Imatrix aset_M(node A, Imatrix M)
{
    node ptc, ptr, ptp, pt0;
    int j, r = 1;

    M = Imatrix_resize(M, aset_npts(A) - aset_r(A), aset_dim(A) - 1);
    ptr = St(A);
    while ((ptc = aset_next_cfg(&ptr)) != 0) {
      pt0 = aset_next_pnt(&ptc);
      while ((ptp = aset_next_pnt(&ptc)) != 0) {
          for (j = 1; j <= aset_dim(A) - 1; j++)
            *IMref(M, r, j) = aset_pnt_get(ptp, j) - aset_pnt_get(pt0, j);
          r++;
      }
    }
    return M;
}

/* end Aset.c */

/**************************************************************************/
/********************** implementations from Types.c **********************/
/**************************************************************************/

static int level=0;
node node_print(node N)
{
  /*DEBUG */
  /*  printf("Entering node_print with node_type = %d.\n", node_get_type(N,LEFT)); */

    if (node_nullp(N) == TRUE)
      return N;
    switch (node_get_type(N, LEFT)) {
    case IDF:
    case STR:
    case ERR:
      printf("%s", (char *) node_get_ptr(N, LEFT));
      break;
    case IMTX:
      Imatrix_print((Imatrix)node_get_ptr(N, LEFT));
      break;
    case DMTX:
      Dmatrix_print((Dmatrix)node_get_ptr(N, LEFT));
      break;
    case PNT:
      pnt_print(N);
      break;
    case PCFG:
      pcfg_print_short(N);
      break;
    case ASET:
      /* aset_print_short(N); */
      aset_print(N);
      break;
    case DBL:
      printf("%f ", node_get_double(N, LEFT));
      break;
    case INT:
      printf("%d ", node_get_int(N, LEFT));
      break;
    case CELL:
        printf("<cly_cell>");
         break;
    case NODE:
        if (level++==0)printf("( ");
        if (node_atomp(Car(N))==FALSE) printf("(");
        node_print((node) Car(N));
        if (node_atomp(Cdr(N))==TRUE&&node_nullp(Cdr(N))==FALSE)
          printf(" . ");
        else printf("  ");
        node_print((node) Cdr(N));
        if ((node_atomp(Cdr(N))==TRUE)||(node_nullp(Cdr(N))==TRUE)){ 
           printf(" )");
        }
        level--;
      break;
    case PROC: printf("%d",(int)node_get_ptr(N,LEFT));
        break;
    default:
      printf("Unknown type %d in Node_Print\n",
             node_get_type(N, LEFT));
      break;
    }
    return N;
}


node atom_int(int val)
{
    node R;
    R = node_new();
    node_set_int(R, val, INT, LEFT);
    return R;
}

node atom_double(double val)
{
    node R;
    R = node_new();
    node_set_double(R, val,DBL, LEFT);
    return R;
}

node atom_id(char *val)
{
    node R;
    R = node_new();
    node_set_ptr(R, mem_strdup(val), IDF, LEFT);
    return R;
}
node atom_proc(node (*prc)(node))
{
    node R;
    R = node_new();
    node_set_ptr(R, (void *)prc, PROC, LEFT);
    return R;
}

node ERRND(char *s)
{ node ans;                                
  ans=node_new();                          
  node_set_ptr(ans,mem_strdup(s),ERR,LEFT);
  return ans;
}


node atom_new(char *val, int tp)
{
    node R;
    R = node_new();
    node_set_ptr(R, val, tp, LEFT);
    return R;
}



void atom_free(node N)
{

    if (node_nullp(N) == TRUE)
      return;
    if (node_atomp(N) != TRUE)
      return;
    switch (node_get_type(N, LEFT)) {
    case ERR:
    case STR:
    case IDF:
      mem_free((char *) node_get_ptr(N, LEFT));
      break;
    case IMTX:
      Imatrix_free((Imatrix)node_get_ptr(N, LEFT));
      break;
    case DMTX:
      Dmatrix_free((Dmatrix)node_get_ptr(N, LEFT));
      break;
    case PNT:
      pnt_free(N);
      break;
    case CELL:
    case ASET:
    case PCFG:
    case INT:
    case NPTR:
      break;
    default:
      printf("Unknown type %d in atom_free\n",
             node_get_type(N, LEFT));
      break;
    }
}

int numericP(int t){
   switch(t){
      case INT: case DBL: case CMPX: case POLY: return TRUE;
                break;
      default: return FALSE;
   }
   return FALSE;
}

/* end Types.c */


/**************************************************************************/
/********************* implementations from globals.c *********************/
/**************************************************************************/

#define IN_GLOBALS_C TRUE

/* - THIS DIDN'T WORK
#define Pel_Out stdout
#define Pel_Err stderr
#define Pel_In  stdin
#define Pel_Log stdout
*/

/*
FILE *Pel_Out;
FILE *Pel_Err;
FILE *Pel_In;
FILE *Pel_Log;
*/

/*
Pel_Out = stdout;
Pel_Err = stderr;
Pel_In  = stdin;
Pel_Log = stdout;
*/
/*
FILE *Pel_Out = stdout;
FILE *Pel_Err = stderr;
FILE *Pel_In  = stdin;
FILE *Pel_Log = stdout;
*/

char *Pel_LogName = 0;
char *FilePrefix =  0;

int Cont_Alg=1;
int Show_Sys=TRUE;
int Show_Xpl=TRUE;

/* end globals.c */

/**************************************************************************/
/********************* implementations from Dtypes.c **********************/
/**************************************************************************/

#define S(i) (DVref(S,i))
void xpnt_unscale(xpnt X,Dvector S){
  int n,i;
  fcomplex ctmp;

  n=(int)xpnt_n(X);
  for(i=1;i<=n;i++){
    ctmp=RCmul(S(i),xpnt_xi(X,i));
    xpnt_xi_set(X,i,ctmp);
  }
}
#undef S
/*
** xpnt_affine
**     rescale an xpnt so that its homogenization coordinate has
**     value 1 (if possible) otherwise its first non-zero coordinate
**     is set to 1.
*/
void xpnt_affine(xpnt X){
 int n,i;
 fcomplex ctmp;
 n=(int)xpnt_n(X);
 if (Cabs(xpnt_h(X))!=0.0){
   for(i=1;i<=n;i++){
      ctmp=Cdiv(xpnt_xi(X,i),xpnt_h(X));
      xpnt_xi_set(X,i,ctmp);
   }
   xpnt_h_set(X,Complex(1.0,0.0));
 }
}

/*
** xpnt_normalize 
**     rescale an xpnt so that it has norm 1 (not including t coord).
*/
void xpnt_normalize(xpnt X){
  int n,i;
  double abs=0.0;
  fcomplex ctmp;
  n=(int)xpnt_n(X);
  
  for(i=1;i<=2*n+2;i++){
    abs+= DVref(X,i)*DVref(X,i);
         
  }
  abs=1/sqrt(abs);
  for(i=1;i<=n;i++){
    ctmp=RCmul(abs,xpnt_xi(X,i));
    xpnt_xi_set(X,i,ctmp);
  }
}

/***********************************************************************/
/***************** implementation code from Extremal.h *****************/
/***********************************************************************/

/*
** global storage for Lin prog solver
*/
int LP_M=0, LP_N=0;
Dmatrix LP_A=0, LP_B=0, LP_C=0, LP_X=0,LP_Q=0, LP_R=0, LP_T1=0, LP_T2=0;
Ivector LP_basis=0, LP_nonbasis=0;
extern double RS_zt;

#define X(i) (DVref(LP_X,i))
#define A(i,j) (DMref(LP_A,i,j))
#define B(i) (DVref(LP_B,i))
#define DtypesC(i) (DVref(LP_C,i))


/*
**  int  pcfg_vertex(node pnt,node PC)
** 
** 
**   return TRUE if pnt is a vertex of PC, 
**          FALSE otherwise (also if pnt not in PC)
**  
**   Notes) a)should check its arguments.
**          b)ratio,LP,and basis should be pointers to matrices
**            of the appropriate type, (they can be 0 to begin
**            with). After completion they will point to the space
**            allocated for the linear program, this allows vertex
**            to be called a number of times without having to 
**            constantly reallocate blocks of space.
**          c)ignores lifting values.
**           
*/

int pcfg_vertex(node pnt,
            node PC)
{
    int i, j, N, D;
    node ptr;

    N = pcfg_npts(PC)-1;
    D = pcfg_dim(PC)-1; /* ignore last coord ... treat it as a lifting value*/

  /* 
  ** initialize basis,nonbasis,X and DtypesC
  */ 
    for (i = 1; i <= N; i++){
      *IVref(LP_nonbasis,i)=i;
      X(i)=0.0;   
      DtypesC(i)=0.0;   
    }
   for (i=1; i<= D+1; i++){
      *IVref(LP_basis,i)=N+i;
      X(N+i)=1.0;   
      DtypesC(N+i)=1.0;   
    }

/* load matrix for linear program */
    ptr = PC;
    for (i = 1; i <= N; i++) {
      ptr = Cdr(ptr);
      if (Car(ptr)==pnt) ptr=Cdr(ptr);
      /* positive and negative coeficients for points */
        for (j = 1; j <= D; j++) {
           A( j, i) = *IVref(pnt_coords(Car(ptr)), j);
        }
        A( D+1,i)=1.0;
    }
    for (i=N+1;i<=D+N+1;i++){
      for (j=1;j<=D+1;j++) A(j,i)=0.0;
    }
    A(D+1,N+1)=1.0;
    B(D+1)=1.0;
    for(i=1;i<=D;i++){
      A(i,N+1+i)=1.0;
      B(i)=*IVref(pnt_coords(pnt),i);
      X(N+1+i)=B(i);
    }
 Rsimp(LP_A,LP_B,LP_C,LP_X,LP_basis,LP_nonbasis,LP_R,LP_Q,LP_T1,LP_T2); 
 for(i=N+1;i<=N+1+D;i++) if (X(i)>RS_zt) return TRUE;
 return FALSE;
}
#undef X
#undef A
#undef B
#undef DtypesC

/*
** pcfg_extremal(node PC)  
**         Input: A point configuration PC.
**         Outut: Same point configuration PC.
**  Side Effects: All non-extremal points are removed from PC.
** 
**      Method: Calls pcfg_vertex on each point to see if it
**              is a vertex and removes those that arent.
**             
**       note:
**       The main loop maintains the following invariant conditions
**                    Car(ptr_nxt) is point being checked
**                    Cdr(ptr_old) = ptr_nxt
*/
node pcfg_extremal(node PC)
{
  node ptr_old, ptr_nxt;
  int N,D;
  
/* allocate space for the linear program and initialize */
    N = pcfg_npts(PC)-1;
    D = pcfg_dim(PC)-1; /* ignore last coord ... treat it as a lifting value*/
    LP_M = D+1;
    LP_N = N+D+1;
    LP_A =  Dmatrix_resize(LP_A, LP_M, LP_N);
    LP_B =  Dmatrix_resize(LP_B, 1, LP_M);
    LP_C =  Dmatrix_resize(LP_C, 1, LP_N);
    LP_X =  Dmatrix_resize(LP_X, 1, LP_N);
    LP_Q =  Dmatrix_resize(LP_Q, LP_M, LP_M);
    LP_R =  Dmatrix_resize(LP_R, LP_M, LP_M);
    LP_T1 = Dmatrix_resize(LP_T1,1, LP_M);
    LP_T2 = Dmatrix_resize(LP_T2,1, LP_M);
    LP_basis    = Imatrix_resize(LP_basis,1,  LP_M);
    LP_nonbasis = Imatrix_resize(LP_nonbasis,1, LP_N-LP_M);


  ptr_nxt=Cdr(PC);
  ptr_old=PC;        
  while(ptr_nxt != 0) {
    if (pcfg_vertex(Car(ptr_nxt), PC) != TRUE){
      pcfg_remove(Car(ptr_nxt),PC,ptr_old); 
    }
    else ptr_old= ptr_nxt; 
    ptr_nxt = Cdr(ptr_nxt);
  }

    Dmatrix_free(LP_A); LP_A=0;
    Dmatrix_free(LP_B); LP_B=0;
    Dmatrix_free(LP_C); LP_C=0;
    Dmatrix_free(LP_X); LP_X=0;
    Dmatrix_free(LP_Q); LP_Q=0;
    Dmatrix_free(LP_R); LP_R=0;
    Dmatrix_free(LP_T1); LP_T1=0;
    Dmatrix_free(LP_T2); LP_T2=0;
    Imatrix_free(LP_basis);    LP_basis=0;
    Imatrix_free(LP_nonbasis); LP_nonbasis=0;
  return PC;
}

node aset_extremal(node A)
{
    node ptr, ptc;
    LOCS(1);
    PUSH_LOC(A);
    ptr = aset_start_cfg(A);
    while ((ptc = aset_next_cfg(&ptr)) != 0) {
        pcfg_extremal(ptc);
    }
    POP_LOCS();
    return A;
}

/***********************************************************************/
/****************** implementation code from RSimp.c *******************/
/***********************************************************************/
                  
/*
** Rsimp   revised simplex method (Using Bland's rule) 
**        and a qr factorization to solve the linear equations
**
**      Adapted from algorithms presented in 
**             Linear Approximations and Extensions                  
**             (theory and algorithms)
**             Fang & Puthenpura
**             Prentice Hall, Engelwood Cliffs NJ (1993)
**      and 
**            Linear Programming
**            Chvatal 
**            Freeman and Company, New York, 1983
** 
**      (developed first in Octave, many thanks to the author)
** 
**
**  Solve the problem 
**       minimize C'x, 
**       subject to A*x=b,  x>=0
**       for x,c,b n-vectors, and A an m,n matrix with full row rank 
** 
** Assumptions:
**    A mxn matrix with full row rank.
**    b an m matrix. 
**    c an n-vector.
**    x an n-vector holding a basic feasible solution, 
**    basis m-vector holding indices of the basic variables in x
**    nonbasis n-m vector holding the indices not appearing in x.
** 
**  Returns: 
**      FAIL if algorithm doesn't terminate.
**      UNBD if problem is unbounded
**      OPT  if optimum found
**  efects:
**    A,b,c unchanged.
**    x basis, nonbasis, hold info describing last basic feasible 
**                       solution.
**    Q,R hold qrdecomp of last basis matrix.
**    t1,t2 undefined.
**
**
*/ 

/* #include "pelutils.h" */

 #define verbose 0
 #define OPT 0
 #define UNBD 1
 #define FAIL -1
 #define zero_tol RS_zt
 #define nonbasis(j) (*IVref(nonbasis,j))
 #define basis(j) (*IVref(basis,j))
 #define AN(i,j) (DMref(A,i,nonbasis(j)))
 #define AB(i,j) (DMref(A,i,basis(j)))
 #define CB(i) (DVref(c,basis(i)))
 #define CN(i) (DVref(c,nonbasis(i)))
 #define XB(i) (DVref(x,basis(i)))
 #define XN(i) (DVref(x,nonbasis(i)))
 #define Y(i)  (DVref(t1,i))
 #define W(i)  (DVref(t2,i))
 #define DtypesD(i)  (DVref(t2,i))
 #define DtypesR(i,j) (DMref(DtypesR,i,j))
 #define Q(i,j) (DMref(Q,i,j))

double RS_zt=0.0000001;

int Rsimp(Dmatrix A, Dvector b, Dvector c,
      Dvector x,Ivector basis,Ivector nonbasis,
      Dmatrix DtypesR, Dmatrix Q, Dvector t1, Dvector t2){
 int m,n,i,j,k,l,q,qv;
 int max_steps=20;
 double r,a,at;

 m=(int)DMrows(A);
 n=(int)DMcols(A);
 max_steps=4*n;

 /*
 **  Dimension assumptions:
 **     A(m,n) x(n) c(n) b(m) Q(mxm) DtypesR(m,m) t1(m) t2(m) 
 **     basis(m) nonbasis(n-m)   
 */   

 for(k=0; k<=max_steps;k++){
   /*
   ** Step 0) load new basis matrix and factor it
   */
    for(i=1;i<=m;i++){
      for(j=1;j<=m;j++){ 
        DtypesR(i,j)=AB(i,j);
      }
    }
    Dmatrix_GQR(Q,DtypesR);

   /*
   ** Step 1) solving system  B'*w=c(basis)
   **      a) forward solve DtypesR'*y=c(basis)
   */
    for(i=1;i<=m;i++){
      Y(i)=0.0;
      for(j=1;j<=i-1;j++){
        Y(i)+=DtypesR(j,i)*Y(j);
      }
      Y(i)=(CB(i)-Y(i))/DtypesR(i,i);
    }
   /*  
   **      b) find w=Q*y
   **         note: B'*w=(Q*DtypesR)'*Q*y= DtypesR'*(Q'*Q)*y=DtypesR'*y=c(basis)
   */
    for(i=1;i<=m;i++){
      W(i)=0.0;
      for(j=1;j<=m;j++){
        W(i)+=Q(i,j)*Y(j);
      }
    }

   /*
   ** Step 2)find entering variable, 
   ** (use lexicographically first variable with negative reduced cost)
   */
    q=n+1;
    for(i=1;i<=n-m;i++){
      /* calculate reduced cost */
      r=CN(i);
      for(j=1;j<=m;j++){
        r-=W(j)*AN(j,i);      
      }
      if (r<-zero_tol){
        if (q==n+1 || nonbasis(i)<nonbasis(q)){
          q=i;                 
        }
      }
    }
   /*
   ** if ratios were all nonnegative current solution is optimal
   */
    if (q==n+1){           
      if (verbose>0) printf("optimal solution found in %d iterations\n",k);
      return OPT;
    }
   /* 
   ** Step 3)Calculate translation direction for q entering
   **        by solving system  B*d=-A(:,nonbasis(q));
   **   a) let y=-Q'*A(:,nonbasis(q));
   */
    for(i=1;i<=m;i++){
      Y(i)=0.0;
      for(j=1;j<=m;j++){
        Y(i)-=Q(j,i)*AN(j,q);
      }
    }
  
   /*
   **  b) back solve Rd=y  (d=R\y)
   **     note B*d= Q*DtypesR*d=Q*y=Q*-Q'*A(:nonbasis(q))=-A(:,nonbasis(q)) 
   */
    for(i=m;i>=1;i--){
      DtypesD(i)=0.0;
      for(j=m;j>=i+1;j--){
        DtypesD(i)+=DtypesR(i,j)*DtypesD(j);
      }
      DtypesD(i)=(Y(i)-DtypesD(i))/DtypesR(i,i);
    }
   /*
   ** Step 4 Choose leaving variable 
   **     (first variable to become negative, by moving in direction DtypesD)
   **     (if none become negative, then objective function unbounded)
   */
    a=0;
    l=0;
    for(i=1;i<=m;i++){
      if (DtypesD(i)<-zero_tol){
         at=-1*XB(i)/DtypesD(i);
         if (l==0 || at<a){
           a=at;
           l=i;
         }
      }
    }
    if (l==0){
      if (verbose>0){
         printf("Objective function Unbounded (%d iterations)\n",k);
      }
      return UNBD;
    }
   /*
   ** Step 5) Update solution and basis data
   */
    XN(q)=a;        
    for(j=1;j<=m;j++){
      XB(j)+=a*DtypesD(j);
    }
    XB(l)=0.0;       /* enforce strict zeroness of nonbasis variables */
    qv=nonbasis(q);
    nonbasis(q)=basis(l);
    basis(l)=qv;
  }
  if (verbose>=0){ 
      printf("Simplex Algorithm did not Terminate in %d iterations\n",k);
  }
  return FAIL;
}

/***********************************************************************/
/******************* implementation code from MSD.c ********************/
/***********************************************************************/

/*
**    copyright (c) 1995  Birk Huber
*/

/*
** MSD.c     
**        Prune and search based computation of cells of mixed subdivisions
*/

/* #include "pelutils.h" -Seems to be redundant */

/* node atom_new(); */
/* node set_up_FaceLists(); */

#define  msd_out  stdout /* was Pel_Out */

/*
*********************************************************************
** Aset_I: internal representation for Asets 
** logical layout:
**   (indices start at one for all vector typesin Aset_I)
**    A_r(A)      ---- number of point configs
**    A_n(A)      ---- dimension (without the lifting value).
**    A_m(A)      ---- total number of points in Aset
**    A_npts(A,i) ---- number of points in ith polytope 
**    A_pt(A,i,j) ---- return jth point of ith polytope
**    A_pt_coords(A,i,j) --- return coordinate vector of jth pnt of ith ptope
**    A_pt_coord(A,i,j,k) --- return kth coord of jth pnt of ith ptope
**
**  physical layout:
**  r,n,m are stored seperatly as integers.
**  all points are stored in one C-vector, with indexing information 
**  stored in one long integer vector brocken up logically as
**  |m_1,....,m_r|s1,....,sr|
*/
#define A_r(A) (A->r)
#define A_n(A) (A->n)
#define A_m(A) (A->m)
#define A_npts(A,r) (((A->store)[(r)-1]))
#define A_ptst(A,r) (((A->store)[A_r(A)+(r)-1]))
#define A_pt(A,r,i) (((A->pts)[A_ptst(A,r)+(i)-1]))
#define A_pt_flat(A,i) (((A->pts)[(i)-1]))
#define A_pt_coords(A,i,j) (pnt_coords(A_pt(A,i,j)))
#define A_pt_coord(A,i,j,k) (*IVref(A_pt_coords(A,i,j),k))

/* NOW IN pelutils.h
typedef struct Aset_Itag{
     int r;
     int n;
     int m;
     node *pts;   */ /*a vector of points */ /*
     int  *store;
     }*Aset_I; */

/*
** Internalize_Aset 
**    Input: an Aset S.
**   Output: an internal (Aset_I) representation of S.
*/
Aset_I internalize_aset(aset S){
  Aset_I A;
  int r=0,pt=0;
  node ptr,ptc,ptp;

  if((A=(Aset_I)malloc(sizeof(struct Aset_Itag)))==0)
      bad_error("malloc failed in internalize_aset");
  A_r(A)=aset_r(S);
  A_n(A)=aset_dim(S)-1;
  A_m(A)=aset_npts(S);
  if((A->store=(int *)malloc((2*A_r(A))*sizeof(int)))==0)
     bad_error("malloc failed in internalize_aset");
  if((A->pts=(node *)malloc(A_m(A)*sizeof(node)))==0)
     bad_error("malloc failed in internalize_aset");

  ptr = aset_start_cfg(S);
  while ((ptc = aset_next_cfg(&ptr)) != 0) {
        A_npts(A,++r)=0;
        A_ptst(A,r)=pt;
        while ((ptp = aset_next_pnt(&ptc)) != 0) {
            A_npts(A,r)++;
            A_pt_flat(A,++pt)=ptp;
        }
    }
  return A;
}
/* 
** Aset_I_free
**     free all storage allocated to an Aset_I.
*/
void Aset_I_free(Aset_I A){
 if (A!=0){
    free((char *)A->store);
    free((char *)A->pts);
    free((char *)A);
 }
}
/*
** Aset_I_print
**    Display internal aset (for debugging purposes)
*/
void Aset_I_print(Aset_I A){
 int i,r;
#ifdef LOG_PRINT   
fprintf(msd_out," R=%d, N=%d, M=%d \n",A_r(A),A_n(A),A_m(A));
   fprintf(msd_out," tops = <")
#endif
;
   for(i=1;i<=A_r(A);i++){
 #ifdef LOG_PRINT    
    fprintf(msd_out,"%d ",A_npts(A,i))
#endif
;
         if (i<A_r(A)) 
#ifdef LOG_PRINT
fprintf(msd_out,",")
#endif
;
   }
 #ifdef LOG_PRINT
  fprintf(msd_out,">\n");
   fprintf(msd_out," sts = <")
#endif
;
      for(i=1;i<=A_r(A);i++){       
#ifdef LOG_PRINT
         fprintf(msd_out,"%d ",A_ptst(A,i))
#endif
;
         if (i<A_r(A)) 
#ifdef LOG_PRINT
fprintf(msd_out,",")
#endif
;
   }
#ifdef LOG_PRINT
   fprintf(msd_out,">\n")
#endif
;
   for(r=1;r<=A_r(A);r++){
#ifdef LOG_PRINT
    fprintf(msd_out,"configuration %d :\n",r)
#endif
;
    for(i=1;i<=A_npts(A,r);i++) {
      pnt_print(A_pt(A,r,i));
    }
   }
} 
/*
*****************************************************************
** Internal representation for mixed cells
** logical layout:
**   C_r(C)          number of point configs
**   C_type(C,i)        dimension of ith peice, 
**                      i from 1 to r.
**   C_idx(C,i,j)     index of jth point of ith peice. 
**                      j from 0 to T[i].      
** 
** physical layout:
**     store[0] ... store[r-1]     hold T.
**     store[r] ... store[2r-1]    hold indices to first rows of Pt
**     store[2r] ... store[3r+n-1] hold Pt
**
*/

/* NOW IN pelutils.h
typedef struct Cell_Itag {
  int r;
  int *store;
}*Cell_I;
*/

#define C_r(C) (((C)->r))
#define C_type(C,i) (((C)->store)[(i)-1])
#define C_st(C,i) (((C)->store)[C_r(C)+(i)-1])
#define C_idx(C,r,i) (((C)->store)[C_st(C,r)+(i)])
#define C_pt(A,C,r,i) (A_pt(A,r,C_idx(C,r,i)))
#define C_pt_coords(A,C,r,i) (pnt_coords(C_pt(A,C,r,i)))
#define C_pt_coord(A,C,r,i,j) (*IVref(C_pt_coords(A,C,r,i),j))

/*
** initialize_cell
**   Input:  an Aset A, and a Type vector T.
**   Output: A blank cell to hold fine mixed type T cells of A
*/

Cell_I initialize_cell(Aset_I A, Ivector T){
 Cell_I C;
 int i,j;
 if ((C=(Cell_I)malloc(sizeof(struct Cell_Itag)))==0) 
     bad_error("malloc failed in initialize_cell");
 C_r(C)=A_r(A);
 if((C->store = (int *)malloc((3*C_r(C)+A_n(A))*sizeof(int)))==0) 
      bad_error("malloc failed in initialize_cell");
 j=2*C_r(C);             /* location of C_Pt(C,1,0) */
 for(i=1;i<=C_r(C);i++){
  C_type(C,i)=*IVref(T,i);  /* copy T into C_type feild */
  C_st(C,i)=j;           /* set location of C_Pt(C,i,0)*/
  j+=C_type(C,i)+1;
 }
/* initialize C_Pt to all zeros (asthetic value only)*/
 for(i=1;i<=C_r(C);i++) for(j=0;j<=C_type(C,i);j++) C_idx(C,i,j)=0;

 return C;
}

/*
** Cell_I_free
**     free all storage allocated to a Cell_I.
*/
void Cell_I_free(Cell_I C){
 if (C!=0){
    free((char *)C->store);
    free((char *)C);
 }
}

/*
** Cell_I_print
**    display a cells points
*/
void Cell_I_print(Cell_I C){
 int i,j;
#ifdef LOG_PRINT
 fprintf(msd_out,"[")
#endif
;
 for(i=1;i<=C_r(C);i++){
   for(j=0;j<=C_type(C,i);j++){
#ifdef LOG_PRINT
      fprintf(msd_out,"%d",C_idx(C,i,j))
#endif
;
      if (j<C_type(C,i)) 
#ifdef LOG_PRINT
fprintf(msd_out,", ")
#endif
;
   }
   if (i<C_r(C)) 
#ifdef LOG_PRINT
fprintf(msd_out,": ")
#endif
;
 }
#ifdef LOG_PRINT
 fprintf(msd_out,"]")
#endif
;
}
/*
*******************************************************************
** Face lists  list of faces of an Aset
**
** Logically: A collection of r lists of Ivectors, the Ivectors in the
**            ith list contain the indices of points in the various 
**            faces of type i in the lower hull of A_i
**            List_Start(i) -- a pointer to the node holding a pointer to the
**                             start of the ith list.
**            List_Ptr(i)   -- a pointer used in iterating through the ith
**                             list.
** Physically: The Variable List_Store holds a pointer to the list of
**             of starting points of the lists (ugh):
**    _______      ______      ______
**   | 0 | | |<---|- | | |<---|- | | |<-----List_Store
**    -----+-      ----+-      ----+-
**         |           |           |
**     List_1      List_2       List_R
**   
**                    ________      ________
**   where List_i ==  | | |  -+---->| | |  -+---> ,,,,,,,
**                    --+-----      --+-----
**                     Face_1       Face_2
*/
 node List_Store=0;
 int List_R=0;
 node *List_Start;
 node *List_Ptrs;
 #define LStart(i)  (List_Start[(i)-1])
 #define LPtr(i)    (List_Ptrs[(i)-1])

/*
** init_Face_List_storeage:
**    Input: A number R of lists to create.
**    Output: FALSE if error occurs, TRUE otherwise.
**    Sets List_R to DtypesR.
**    Creates a chain of DtypesR nodes starting at List_Store, (with nodes linked 
**       through their Left pointers, violating the usual convention for lists)
**    The List_Start, and List_Ptr, arrays are created and initialized so that
**       the ith entrees point to the ith nodes on the list List_Store.
**    The Right Pointer of each of these nodes should be interpereted as a 
**       pointer to the start of the corresponding list, initialized to 0.
*/
int init_Face_List_storeage(int r){
 node tmp;
 int i;
 LOCS(1);
 PUSH_LOC(List_Store);
 List_R=r;
 List_Start=(node *)mem_malloc(r*sizeof(node));
 List_Ptrs=(node *)mem_malloc(r*sizeof(node));
 for(i=1;i<=r;i++){
    tmp=node_new();
    node_set_ptr(tmp,(void *)List_Store,NODE,LEFT);
    node_set_ptr(tmp,(void *)0         ,NODE,RIGHT);
    List_Store=tmp;
    LStart(i)=tmp;
    LPtr(i)=tmp;
 }
 POP_LOCS();
 return TRUE;
}
/*
** free_Face_list -- free data storage for the List_Start and List_Ptrs arrays
**                   reset List_Store and List_R to 0
*/
void free_Face_list(){
 mem_free((char *)List_Start);
 mem_free((char *)List_Ptrs);
 List_Store=0;
 List_R=0;
}

/*
** display all faces on face list
*/
void print_Face_list(){
 int i;
 node ptr;
 for(i=1;i<=List_R;i++){
   ptr=LStart(i);
   while((ptr=(node)node_get_ptr(ptr,RIGHT))!=0){
       Imatrix_print((Imatrix)node_get_ptr(ptr,LEFT));
#ifdef LOG_PRINT     
  fprintf(msd_out,"\n")
#endif
;
   }
 }
}

/*
** Push_Face:
*/
void Push_Face(Imatrix F,int r){
    node tmp,ptr;
    tmp=node_new();
    ptr=(node)node_get_ptr(LPtr(r),RIGHT);
    node_set_ptr(tmp,(void *)F,IMTX,LEFT);
    node_set_ptr(tmp,(void *)ptr,NODE,RIGHT);
    node_set_ptr(LPtr(r),(void *)tmp,NODE,RIGHT);
}

void zeroth_face(int j){
 LPtr(j)=LStart(j);
}

int next_face(Cell_I *C,Aset_I A, int r){
 int i;
 Imatrix F;
 /* advance pointer */
 if (LPtr(r)==0) return 0;
 LPtr(r)=(node)node_get_ptr(LPtr(r),RIGHT);
 if (LPtr(r)==0) return 0;
 /* copy face */
 F=(Imatrix)node_get_ptr(LPtr(r),LEFT);
 for(i=0;i<=C_type(*C,r);i++) C_idx(*C,r,i)=*IVref0(F,i);
 return 1;
}

/*
** Intermediate testing -- handles testing for all incomplete cells
*/

static int MSD_LP_M=0, MSD_LP_N=0;
static Dmatrix MSD_LP_A=0, MSD_LP_B=0, MSD_LP_C=0, MSD_LP_X=0;
static Dmatrix MSD_LP_Q=0, MSD_LP_R=0, MSD_LP_T1=0, MSD_LP_T2=0;
static Ivector MSD_LP_basis=0, MSD_LP_nonbasis=0;

/*
** set_up_LP
**   Input: An Aset_I A and a Cell_I C 
**   Output: LP,Ratio, and Basis have enough space reserved for them
**           to allow them to hold any of the linear programs 
**           required by MSD for A and C with only reseting bounds.
*/
void set_up_LP(Aset_I A, Cell_I C){
  MSD_LP_M=A_n(A)+A_r(A);
  MSD_LP_N=A_m(A);
  MSD_LP_A =  Dmatrix_resize(MSD_LP_A, MSD_LP_M, MSD_LP_N);
  MSD_LP_B =  Dmatrix_resize(MSD_LP_B, 1, MSD_LP_M);
  MSD_LP_C =  Dmatrix_resize(MSD_LP_C, 1, MSD_LP_N);
  MSD_LP_X =  Dmatrix_resize(MSD_LP_X, 1, MSD_LP_N);
  MSD_LP_Q =  Dmatrix_resize(MSD_LP_Q, MSD_LP_M, MSD_LP_M);
  MSD_LP_R =  Dmatrix_resize(MSD_LP_R, MSD_LP_M, MSD_LP_M);
  MSD_LP_T1 = Dmatrix_resize(MSD_LP_T1,1, MSD_LP_M);
  MSD_LP_T2 = Dmatrix_resize(MSD_LP_T2,1, MSD_LP_M);
  MSD_LP_basis    = Imatrix_resize(MSD_LP_basis,1,  MSD_LP_M);
  MSD_LP_nonbasis = Imatrix_resize(MSD_LP_nonbasis,1, MSD_LP_N-MSD_LP_M);
}

void free_LP(){
  MSD_LP_M=0; MSD_LP_N=0;
  Dmatrix_free(MSD_LP_A); MSD_LP_A=0;
  Dmatrix_free(MSD_LP_B); MSD_LP_B=0;
  Dmatrix_free(MSD_LP_C); MSD_LP_C=0;
  Dmatrix_free(MSD_LP_X); MSD_LP_X=0;
  Dmatrix_free(MSD_LP_Q); MSD_LP_Q=0;
  Dmatrix_free(MSD_LP_R); MSD_LP_R=0;
  Dmatrix_free(MSD_LP_T1); MSD_LP_T1=0;
  Dmatrix_free(MSD_LP_T2); MSD_LP_T2=0;
  Imatrix_free(MSD_LP_basis);    MSD_LP_basis=0;
  Imatrix_free(MSD_LP_nonbasis); MSD_LP_nonbasis=0;
}

/*
** Load_LP
**
*/

#define MSD_A(i,j)      (DMref(MSD_LP_A,i,j))
#define MSD_B(i)        (DVref(MSD_LP_B,i))
#define MSD_C(i)        (DVref(MSD_LP_C,i))
#define MSD_X(i)        (DVref(MSD_LP_X,i))
#define MSDbasis(i)    (*IVref(MSD_LP_basis,i))
#define MSDnonbasis(i) (*IVref(MSD_LP_nonbasis,i))

void Load_LP(int st,int tp,Cell_I C, Aset_I A){
   int cidx,bidx;
   int i,j,k;
   int nrows,ncols;
   int n;          /* dimenstion of points (without lifting) */
   int tm=0;       /* total number of points in configs ts..tp  */
   int tr=1+tp-st; /* number of distinct configs to consider*/
   int tc=0;
   double tmp;
   Imatrix coords;

   /* calculate new matrix dims */
   n=A_n(A);
   for(i=st;i<=tp;i++){
      tm+=A_npts(A,i);
   }
   nrows=n+tr;
   ncols=tm+n;
   MSD_LP_A=Dmatrix_resize(MSD_LP_A,nrows,ncols);
   MSD_LP_B=Dmatrix_resize(MSD_LP_B,1,nrows);
   MSD_LP_C=Dmatrix_resize(MSD_LP_C,1,ncols);
   MSD_LP_X=Dmatrix_resize(MSD_LP_X,1,ncols);
   MSD_LP_Q=Dmatrix_resize(MSD_LP_Q,nrows,nrows);
   MSD_LP_R=Dmatrix_resize(MSD_LP_R,nrows,ncols);
   MSD_LP_T1=Dmatrix_resize(MSD_LP_T1,1,nrows);
   MSD_LP_T2=Dmatrix_resize(MSD_LP_T2,1,nrows);
   MSD_LP_basis=Imatrix_resize(MSD_LP_basis,1,nrows);
   MSD_LP_nonbasis=Imatrix_resize(MSD_LP_nonbasis,1,ncols-nrows);

   for(j=1;j<=ncols;j++){
        MSD_X(j)=0.0;
        for(i=1;i<=nrows;i++) MSD_A(i,j)=0;
   }
   for(i=1;i<=nrows;i++) MSD_B(i)=0; 

   /* 
   ** put n independent collumns at start of lp
   */
   for(j=1;j<=n;j++){
     for(i=1;i<=nrows;i++) MSD_A(i,j)=rand_double(0,10);
     MSD_C(j)=1500.0;
   }

   /* 
   ** 
   */
   cidx=n;
   bidx=0;
   for(i=st;i<=tp;i++){
      for(j=1;j<=A_npts(A,i);j++){
         MSD_A(1+i-st,cidx+j)=1.0;
         coords=pnt_coords(A_pt(A,i,j));
         for(k=1;k<=n;k++){
            MSD_A(tr+k,cidx+j)=*IVref(coords,k);
         }
         MSD_C(cidx+j)=*IVref(coords,n+1); 
      }
      MSD_B(1+i-st)=1.0;
      tmp=1.0/(C_type(C,i)+1);
      for(j=0;j<=C_type(C,i);j++){
         MSDbasis(++bidx)=cidx+C_idx(C,i,j);
         MSD_X(cidx+C_idx(C,i,j))=tmp;
         for(k=1;k<=n;k++){
           MSD_B(tr+k)+=tmp*MSD_A(tr+k,cidx+C_idx(C,i,j));
         }
      }
   cidx+=A_npts(A,i);
   }
   for(j=1;j<=n;j++){
      if (j<=nrows-bidx) MSDbasis(bidx+j)=j;
      else MSDnonbasis(j-nrows+bidx)=j;
   }
   bidx=bidx-nrows+n;
   for (i=n-tc+1;i<=ncols;i++) if(MSD_X(i)==0.0)MSDnonbasis(++bidx)=i;
}
#undef MSD_A
#undef MSD_B
#undef MSD_C
#undef MSD_X
#undef MSDbasis
#undef MSDnonbasis

/*
** Final testing -- Verification of complete cells
*/
Imatrix M;
Imatrix U; 
Imatrix Norm;
int vol;
void set_up_Final(int n){
 M=Imatrix_new(n,n+1);
 U=Imatrix_new(n,n);
 Norm=Imatrix_new(1,n+1);
 vol=0;
}

void free_Final(){
 Imatrix_free(M);
 Imatrix_free(U);
 Imatrix_free(Norm);
 M=0;Norm=0;U=0;vol=0;
}


int Final_Check(Aset_I A, Cell_I C){
 int i,j,k;
 int row=0,s0,s1;
 
 /* load matrix */
 for(i=1;i<=A_r(A);i++){
    for(j=1;j<=C_type(C,i);j++){
       row++;
       for(k=1;k<=A_n(A)+1;k++){
         *IMref(M,row,k)=C_pt_coord(A,C,i,j,k)-C_pt_coord(A,C,i,0,k);
       }
    }
 }
 /* factor matrix */
 Imatrix_hermite(M,U);

 /* calculate volume */
 vol=1;
 for(i=1;i<=A_n(A);i++) (vol)*=*IMref(M,i,i);

 if (vol==0) return FALSE;

 /* calculate normal and fix direction */
 Imatrix_backsolve(M,Norm);
 if (*IVref(Norm,A_n(A)+1)<0){
  for(i=1;i<=A_n(A)+1;i++) *IVref(Norm,i)*=-1;
 } 
 
 for(i=1;i<=A_r(A);i++){
    /* find offset of first point of cell */
    s0=0;
    for(k=1;k<=A_n(A)+1;k++){
      s0+=*IVref(Norm,k)*C_pt_coord(A,C,i,0,k);
    }
    /*check that remaining points have same norm */
    for(j=1;j<=C_type(C,i);j++){
       s1=0;  
       for(k=1;k<=A_n(A)+1;k++){
           s1+=*IVref(Norm,k)*C_pt_coord(A,C,i,j,k);
       }
       if (s1!=s0){warning("bad normall");
                    return FALSE;}
    }
    /* check that no points have lower offset */
    for(j=1;j<=A_npts(A,i);j++){
       s1=0;                
       for(k=1;k<=A_n(A)+1;k++){
         s1+=(*IVref(Norm,k))*A_pt_coord(A,i,j,k);
       }
       if (s1<s0) return FALSE;
    }
 }

return TRUE;
}


/*
*******************************************************
**  Algorithm: Finding All type t mixed cells
**
**  Input A_1,...,A_r  a lifted Aset
**        t_1,...,t_r  type of cell desired
**
**  Output S a list of all type t cells of our subdivision.
**
**  Begin
**     S={}
**     C_1 = 0
**     i=1;
**     while(i>0){
**       if ((C_i = next_face(C_i,A_i)!=0) {
**         if (i < r){
**            Add rows to simplex tableau to look for C_i
**            if (new system is feasible) {
**               i++;
**               C_i = 0;
**            }
**         }
**         else {
**           if (C is full dim and C lies on lower hull ){
**              S=S union {C_1,...C_n}
**           }
**         }
**       }
**       else i--
**     }
*/

node MSD(aset Ast, Ivector T){
  Aset_I A;                  /* internal representation of Ast*/
  Cell_I C;                  /* internal rep of target cell */
  int i=1;                   /* how many fields of C speceifed*/
  node Normals=0;            /* output normals list */
  int mv=0;                  /* total mixed volume */
  LOCS(2);
  PUSH_LOC(Normals);
  PUSH_LOC(List_Store);

  /* input verification should check A and T for consistancy*/
  if (T==0) bad_error("type vector must be specified in MSD");
  
  /*
  **  Initialization 
  */
  A=internalize_aset(Ast);
  C=initialize_cell(A,T);
  set_up_LP(A,C);
  set_up_Final(A_n(A));
  set_up_FaceLists(A,C);

  /* 
  ** Iteration
  */
  zeroth_face(1);
  while(i>0){
   if (next_face(&C,A,i)>0){
     if (i<A_r(A)){
       Load_LP(1,i,C,A);
       if (TRUE==IsLower()){
         i++;
         zeroth_face(i);
       }
     }
     else {
       if (Final_Check(A,C)==TRUE){
         Cell_I_print(C); 
#ifdef LOG_PRINT
fprintf(msd_out," vol %d\n",vol)
#endif
;
         mv+=vol;
         Normals=Cons(atom_new((char *)Imatrix_dup(Norm,0),IMTX),
                      Normals);         
       }
     }
   }
   else i--;
  }

/* 
**Clean Up
*/
 Aset_I_free(A);
 Cell_I_free(C);
 free_LP();
 free_Final();
 free_Face_list();
 POP_LOCS();

 #ifdef LOG_PRINT  
 fprintf(msd_out,"Mixed Volume = %d",mv)
#endif
;
 return Normals;
 }

/* 
** Enumerate all k_i Faces 
*/

void zeroth_set(Cell_I *C,int j){
 int i;
 for(i=0;i<=C_type(*C,j);i++) C_idx(*C,j,i)=i+1;
 C_idx(*C,j,C_type(*C,j))--;
}


int next_set(Cell_I *C,Aset_I A, int j){
  int cid,top;
  cid=(top=C_type(*C,j));
  while(cid>=0 && cid<=top){
    if (C_idx(*C,j,cid)<(A_npts(A,j)-top+cid)){
        C_idx(*C,j,cid)++;
        if (cid<top) C_idx(*C,j,cid+1)=C_idx(*C,j,cid);
        cid++;
    }
    else cid--;
  }
  return cid;
}


node set_up_FaceLists(Aset_I A, Cell_I C){
  int i,j,k;     
  Imatrix ML,UL,Itmp;
  LOCS(1);
  PUSH_LOC(List_Store);
  ML=Imatrix_new(A_n(A)+1,A_n(A)+1);
  UL=Imatrix_new(A_n(A)+1,A_n(A)+1);

  init_Face_List_storeage(A_r(A));

  for(i=1;i<=A_r(A);i++){
    ML=Imatrix_submat(ML,A_n(A),C_type(C,i));
    UL=Imatrix_submat(UL,A_n(A),A_n(A));
    zeroth_set(&C,i);
    while(next_set(&C,A,i)>0){
      Load_LP(i,i,C,A); 
      if ( TRUE==IsLower()){
        for(j=1;j<=C_type(C,i);j++){
          for(k=1;k<=A_n(A);k++){
            *IMref(ML,k,j)=C_pt_coord(A,C,i,j,k)-
                           C_pt_coord(A,C,i,0,k);
          }
        }
        Imatrix_hermite(ML,UL);
        k=1; for(j=1;j<=C_type(C,i);j++) k*=*IMref(ML,j,j);
        if (k!=0) {
          Itmp=Imatrix_new(1,C_type(C,i)+1);
          for(j=0;j<C_type(C,i)+1;j++)*IVref0(Itmp,j)=C_idx(C,i,j);
          Push_Face(Itmp,i);
        }
      }  
    }   
  }
  Imatrix_free(ML);
  Imatrix_free(UL);
  POP_LOCS();
  return List_Store;
}



 #define verbose 0
 #define OPT 0
 #define UNBD 1
 #define FAIL -1
 #define MSDzero_tol 0.0000001
 #define MSDnonbasis(j) (*IVref(MSD_LP_nonbasis,j))
 #define MSDbasis(j) (*IVref(MSD_LP_basis,j))
 #define MSD_AN(i,j) (DMref(MSD_LP_A,i,MSDnonbasis(j)))
 #define MSD_AB(i,j) (DMref(MSD_LP_A,i,MSDbasis(j)))
 #define MSD_CB(i) (DVref(MSD_LP_C,MSDbasis(i)))
 #define MSD_CN(i) (DVref(MSD_LP_C,MSDnonbasis(i)))
 #define MSD_XB(i) (DVref(MSD_LP_X,MSDbasis(i)))
 #define MSD_XN(i) (DVref(MSD_LP_X,MSDnonbasis(i)))
 #define MSD_Y(i)  (DVref(MSD_LP_T1,i))
 #define MSD_W(i)  (DVref(MSD_LP_T2,i))
 #define MSD_D(i)  (DVref(MSD_LP_T2,i))
 #define MSD_R(i,j) (DMref(MSD_LP_R,i,j))
 #define MSD_Q(i,j) (DMref(MSD_LP_Q,i,j))

int IsLower(){
 int m,n,i,j,k,l,q,qv;
 int max_steps=20;
 double r,a,at;

 m=(int)DMrows(MSD_LP_A);
 n=(int)DMcols(MSD_LP_A);
 max_steps=4*n;

 /*
 **  Dimension assumptions:
 **     A(m,n) x(n) c(n) b(m) Q(mxm) R(m,m) t1(m) t2(m) 
 **     MSDbasis(m) MSDnonbasis(n-m)   
 */   
/*
   printf("\n A=\n");
   Dmatrix_print(MSD_LP_A);printf("\n B=");fflush(stdout);
   Dmatrix_print(MSD_LP_B);printf("\n C=");fflush(stdout);
   Dmatrix_print(MSD_LP_C);printf("\n X=");fflush(stdout);
   Dmatrix_print(MSD_LP_X);printf("\n Basis=");fflush(stdout);
   Imatrix_print(MSD_LP_basis);printf("\n NonBasis=");fflush(stdout);
   Imatrix_print(MSD_LP_nonbasis);printf("\n");
   printf("\n\n");
*/

 for(k=0; k<=max_steps;k++){
   /*
   ** Step 0) load new basis matrix and factor it
   */
    for(i=1;i<=m;i++){
      for(j=1;j<=m;j++){ 
        MSD_R(i,j)=MSD_AB(i,j);
      }
    }
    Dmatrix_GQR(MSD_LP_Q,MSD_LP_R);

   /*
   ** Step 1) solving system  B'*w=c(basis)
   **      a) forward solve R'*y=c(basis)
   */
    for(i=1;i<=m;i++){
      MSD_Y(i)=0.0;
      for(j=1;j<=i-1;j++){
        MSD_Y(i)+=MSD_R(j,i)*MSD_Y(j);
      }
      if (MSD_R(i,i)!=0) MSD_Y(i)=(MSD_CB(i)-MSD_Y(i))/MSD_R(i,i);
      else {
        printf("singular matrix found\n");
        return TRUE;
      }
    }
   /*  
   **      b) find w=Q*y
   **         note: B'*w=(Q*R)'*Q*y= R'*(Q'*Q)*y=R'*y=c(basis)
   */
    for(i=1;i<=m;i++){
      MSD_W(i)=0.0;
      for(j=1;j<=m;j++){
        MSD_W(i)+=MSD_Q(i,j)*MSD_Y(j);
      }
    }

   /*
   ** Step 2)find entering variable, 
   ** (use lexicographically first variable with negative reduced cost)
   */
    q=n+1;
    for(i=1;i<=n-m;i++){
      /* calculate reduced cost */
      r=MSD_CN(i);
      for(j=1;j<=m;j++){
        r-=MSD_W(j)*MSD_AN(j,i);      
      }
      if (r<-MSDzero_tol){
        if (q==n+1 || MSDnonbasis(i)<MSDnonbasis(q)){
          q=i;                 
        }
      }
    }
   /*
   ** if ratios were all nonnegative current solution is optimal
   */
    if (q==n+1){           
      if (verbose>0){
         printf("optimal solution found in %d iterations\n",k);
      }
      return TRUE;
    }
   /* 
   ** Step 3)Calculate translation direction for q entering
   **        by solving system  B*d=-A(:,MSDnonbasis(q));
   **   a) let y=-Q'*A(:,MSDnonbasis(q));
   */
    for(i=1;i<=m;i++){
      MSD_Y(i)=0.0;
      for(j=1;j<=m;j++){
        MSD_Y(i)-=MSD_Q(j,i)*MSD_AN(j,q);
      }
    }
  
   /*
   **  b) back solve Rd=y  (d=R\y)
   **     note B*d= Q*R*d=Q*y=Q*-Q'*A(:MSDnonbasis(q))=-A(:,MSDnonbasis(q)) 
   */
    for(i=m;i>=1;i--){
      MSD_D(i)=0.0;
      for(j=m;j>=i+1;j--){
        MSD_D(i)+=MSD_R(i,j)*MSD_D(j);
      }
      MSD_D(i)=(MSD_Y(i)-MSD_D(i))/MSD_R(i,i);
    }
   /*
   ** Step 4 Choose leaving variable 
   **     (first variable to become negative, by moving in direction D)
   */
    a=0;
    l=0;
    for(i=1;i<=m;i++){
      if (MSD_D(i)<-MSDzero_tol){
         at=-1*MSD_XB(i)/MSD_D(i);
         if (l==0 || at<a){
           a=at;
           l=i;
         }
      }
    }
    if (l==0){
      if (verbose>0){
         printf("Objective function Unbounded (%d iterations)\n",k);
      }
      return FAIL;
    }
   /*
   ** Step 5) If step is non-degenerate stop. otherwise update basis 
   */
   if (a>=MSDzero_tol){
     if (verbose>0){
       printf("non-degenerate step after %d iterations\n",k);
     }
     return FALSE;
    }
    qv=MSDnonbasis(q);
    MSDnonbasis(q)=MSDbasis(l);
    MSDbasis(l)=qv;
   }
   
  if (verbose>=0){ 
      printf("Simplex Algorithm did not Terminate in %d iterations\n",k);
  }
  return FAIL;
}



/* end Dtypes.c */




Generated by  Doxygen 1.6.0   Back to index