Changeset c7f3b7 in git

Timestamp:

Dec 3, 1999, 12:03:52 AM (24 years ago)

Author:

Moritz Wenk <wenk@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

9f3e9bef10a5b0f7fd4ed74bc83407329b62b4a6

Parents:

2983b3377c3bd9abf59dd57780cedc5d940ca219

Message:

wenk: added new commands: simplex, convhull (mit Greuel abgesprochen)


git-svn-id: file:///usr/local/Singular/svn/trunk@3951 2c84dea3-7e68-4137-9b89-c4e89433aadc

Location:

Singular

Files:

• LIB/solve.lib (modified) (4 diffs)
• iparith.cc (modified) (5 diffs)
• mpr_base.cc (modified) (51 diffs)
• mpr_base.h (modified) (2 diffs)
• mpr_global.h (modified) (2 diffs)
• mpr_inout.cc (modified) (13 diffs)
• mpr_inout.h (modified) (2 diffs)
• mpr_numeric.cc (modified) (22 diffs)
• mpr_numeric.h (modified) (2 diffs)
• tok.h (modified) (3 diffs)

Singular/LIB/solve.lib

@@ -1 +1 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-version="$Id: solve.lib,v 1.13 1999-11-14 21:35:02 wenk Exp $";
+version="$Id: solve.lib,v 1.14 1999-12-02 23:03:52 wenk Exp $";
 info="
 LIBRARY: solve.lib     PROCEDURES TO SOLVE POLYNOMIAL SYSTEMS
@@ -62 +62 @@
         {
             ERROR("Third parameter l must be positive!");
-        }
-    }
+        }    }
     if ( size(#) >= 3 )
     {
@@ -883 +882 @@
 }
 
-
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -946 +944 @@
 }
 
+///////////////////////////////////////////////////////////////////////////////
+
+proc simplexOut(list l)
+"USAGE:   simpelxOut(l); l list
+ASSUME:
+RETURN:
+EXAMPLE: example simplexOut; shows an example
+"
+{
+  int i,j;
+  matrix m= l[1];
+  intvec iposv= l[3];
+  int icase= l[2];
+
+  int cols= ncols(m);
+  int rows= nrows(m);
+
+  int N= l[6];
+
+  if ( -1 == icase )  // objective function is unbound
+  {
+    "objective function is unbound";
+    return;
+  }
+  if ( 1 == icase )  // no solution satisfies the given constraints
+  {
+    "no solution satisfies the given constraints";
+    return;
+  }
+  if ( -2 == icase )  // other error
+  {
+    "an error occurred during simplex computation!";
+    return;
+  }
+
+  for ( i = 1; i <= rows; i++ )
+  {
+    if (i == 1) 
+    {
+      "z = "+string(m[1][1]);
+    }
+    else 
+    {
+      if ( iposv[i-1] <= N+1 ) 
+      {
+        "x"+string(i-1)+" = "+string(m[1][iposv[i-1]]);
+      }
+//        else 
+//        {
+//           "Y"; iposv[i-1]-N+1;
+//        }
+    }
+  } 
+}
+example
+{
+    "EXAMPLE:";echo=2;
+    ring r = (real,10),(x),lp;
+
+    // suppose we have the 
+    
+    matrix sm[5][5]=(  0, 1, 1, 3,-0.5,
+                     740,-1, 0,-2, 0,
+                       0, 0,-2, 0, 7,
+                     0.5, 0,-1, 1,-2,
+                       9,-1,-1,-1,-1);
+
+    // simplex input:
+    // 1: matrix
+    // 2: number of variables
+    // 3: total number of constraints (#4+#5+#6)
+    // 4: number of <= constraints
+    // 5: number of >= constraints
+    // 6: number of == constraints
+
+    list sol= simplex(sm, 4, 4, 2, 1, 1);
+    simplexOut(sol);
+}
+
+///////////////////////////////////////////////////////////////////////////////
 
 // local Variables: ***

Singular/iparith.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: iparith.cc,v 1.193 1999-12-01 13:22:52 Singular Exp $ */
+/* $Id: iparith.cc,v 1.194 1999-12-02 23:03:48 wenk Exp $ */
 
 /*
@@ -138 +138 @@
   { "continue",    0, CONTINUE_CMD ,      CONTINUE_CMD},
   { "contract",    0, CONTRACT_CMD ,      CMD_2},
+  { "convhull",    0, NEWTONPOLY_CMD,     CMD_1},
   { "dbprint",     0, DBPRINT_CMD ,       CMD_M},
   { "def",         0, DEF_CMD ,           ROOT_DECL},
@@ -278 +279 @@
   { "rvar",        0, IS_RINGVAR ,        CMD_1},
   { "setring",     0, SETRING_CMD ,       SETRING_CMD},
+  { "simplex",     0, SIMPLEX_CMD,        CMD_M},
   { "simplify",    0, SIMPLIFY_CMD ,      CMD_2},
   { "size",        0, COUNT_CMD ,         CMD_1},
@@ -3612 +3614 @@
 ,{kWeight,      WEIGHT_CMD,      INTVEC_CMD,     MODUL_CMD }
 ,{jjLOAD1,      LOAD_CMD,        NONE,           STRING_CMD }
+,{loNewtonP,    NEWTONPOLY_CMD,  IDEAL_CMD,      IDEAL_CMD}
 ,{NULL,         0,               0,              0}
 };
@@ -4913 +4916 @@
 ,{jjUNLOAD,    UNLOAD_CMD,      NONE,               -2 }
 #endif /* HAVE_NAMESPACES */
+,{loSimplex,   SIMPLEX_CMD,     LIST_CMD,            6 }
 ,{nuUResSolve, URSOLVE_CMD,     LIST_CMD,            4 }
 ,{NULL,        0,               0,                   0 }

Singular/mpr_base.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: mpr_base.cc,v 1.15 1999-11-23 17:58:47 Singular Exp $ */
+/* $Id: mpr_base.cc,v 1.16 1999-12-02 23:03:49 wenk Exp $ */
 
 /*
@@ -9 +9 @@
  */
 
+#include <math.h>
 #include <limits.h>
-#include <math.h>
 
 #include "mod2.h"
@@ -26 +26 @@
 #include "numbers.h"
 #include "longalg.h"
-#include "tok.h"
 #ifdef HAVE_FACTORY
 #include "clapsing.h"
@@ -43 +42 @@
 #define FreeSizeOf(X,Y) Free(X,sizeof(Y))
 #endif
-
 //<-
 
@@ -70 +68 @@
 //-> sparse resultant matrix
 
-//-> typedefs and structs
-
 /* set of points */
 class pointSet;
 
-/* Linear Program stuff */
-struct linProg
-{
-  mprfloat **LiPM;
-  int *izrov, *iposv;
-  int LiPM_cols,LiPM_rows;
-};
-
-typedef linProg *linProgP;
+
 
 /* sparse resultant matrix class */
@@ -144 +132 @@
   ideal rmat;        // sparse matrix representation
 
-  linProg LP;
+  simplex * LP;      // linear programming stuff
 };
 //<-
 
+//-> typedefs and structs
+poly monomAt( poly p, int i );
 
 typedef unsigned int Coord_t;
@@ -270 +260 @@
 {
 public:
-  convexHull( linProgP _pLP ) : pLP(_pLP) {}
+  convexHull( simplex * _pLP ) : pLP(_pLP) {}
   ~convexHull() {}
 
@@ -277 +267 @@
    * Returns Q[].
    */
-   pointSet ** newtonPolytopes( ideal gls );
+  pointSet ** newtonPolytopesP( const ideal gls );
+  ideal newtonPolytopesI( const ideal gls );
 
 private:
@@ -288 +279 @@
   pointSet **Q;
   int n;
-  linProgP pLP;
+  simplex * pLP;
 };
 //<-
@@ -297 +288 @@
 {
 public:
-  mayanPyramidAlg( linProgP _pLP ) : n(pVariables), pLP(_pLP) {}
+  mayanPyramidAlg( simplex * _pLP ) : n(pVariables), pLP(_pLP) {}
   ~mayanPyramidAlg() {}
 
@@ -342 +333 @@
   Coord_t acoords[MAXVARS+2];
 
-  linProgP pLP;
+  simplex * pLP;
 };
-
-poly monomAt( poly p, int i );
+//<-
 
 #ifndef ASO_GENERATE
-//<-
 
 //-> debug output stuff
@@ -441 +430 @@
   for ( i= 0; i <= max; i++ )
   {
-    points[i]= (onePointP)AllocSizeOf( onePoint );
+    points[i]= (onePointP)Alloc( sizeof(onePoint) );
     points[i]->point= (Coord_t *)Alloc0( (dim+2) * sizeof(Coord_t) );
   }
@@ -454 +443 @@
   {
     Free( (ADDRESS) points[i]->point, fdim * sizeof(Coord_t) );
-    FreeSizeOf( (ADDRESS) points[i], onePoint );
+    Free( (ADDRESS) points[i], sizeof(onePoint) );
   }
   Free( (ADDRESS) points, (max+1) * sizeof(onePointP) );
@@ -476 +465 @@
     for ( i= max+1; i <= max*2; i++ )
     {
-      points[i]= (onePointP)AllocSizeOf( onePoint );
+      points[i]= (onePointP)Alloc( sizeof(struct onePoint) );
       points[i]->point= (Coord_t *)Alloc0( fdim * sizeof(Coord_t) );
     }
@@ -720 +709 @@
   }
 
-#ifdef mprDEBUG_PROT
+#ifdef mprDEBUG_ALL
   Print(" lift vector: ");
   for ( j=1; j < dim; j++ ) Print(" %d ",l[j] );
@@ -754 +743 @@
 bool convexHull::inHull(poly p, poly pointPoly, int m, int site)
 {
-  int i, j, col, icase;
-  int numcons;                // num of constraints
-  int numpts;                // num of pts in defining support
-  int numcols;                // tot number of cols
-
-  numcons = n+1;
-  numpts = m-1;
-  numcols = numpts+1;                // this includes col of cts
-
-  pLP->LiPM[1][1] = +0.0;  pLP->LiPM[1][2] = +1.0;        // optimize (arbitrary) var
-  pLP->LiPM[2][1] = +1.0;  pLP->LiPM[2][2] = -1.0;         // lambda vars sum up to 1
-  for ( j=3; j<=numcols; j++)
-  {
-    pLP->LiPM[1][j] = +0.0; pLP->LiPM[2][j] = -1.0;
+  int i, j, col;
+
+  pLP->m = n+1;
+  pLP->n = m;                // this includes col of cts
+
+  pLP->LiPM[1][1] = +0.0;  
+  pLP->LiPM[1][2] = +1.0;        // optimize (arbitrary) var
+  pLP->LiPM[2][1] = +1.0;  
+  pLP->LiPM[2][2] = -1.0;         // lambda vars sum up to 1
+
+  for ( j=3; j <= pLP->n; j++)
+  {
+    pLP->LiPM[1][j] = +0.0; 
+    pLP->LiPM[2][j] = -1.0;
   }
 
@@ -785 +774 @@
 #ifdef mprDEBUG_ALL
   Print("Matrix of Linear Programming\n");
-  print_mat( pLP->LiPM, numcons+1, numcols);
-#endif
-#if 1
-  if ( numcons + 1 > pLP->LiPM_rows )
-    WerrorS("convexHull::inHull: #rows > #pLP->LiPM_rows!");
-  if ( numcols + 1 > pLP->LiPM_cols )
-    WerrorS("convexHull::inHull: #cols > #pLP->LiPM_cols!");
-#endif
-
-  simplx( pLP->LiPM, numcons, numcols-1, 0, 0, numcons, &icase, pLP->izrov, pLP->iposv);
-
-  return (icase == 0);
+  print_mat( pLP->LiPM, pLP->m+1,pLP->n);
+#endif
+
+  pLP->m3= pLP->m;
+
+  pLP->compute();
+
+  return (pLP->icase == 0);
 }
 
@@ -802 +787 @@
 // ST_SPARSE_VADD: new vertex of convex hull added
 // ST_SPARSE_VREJ: point rejected (-> inside hull)
-pointSet ** convexHull::newtonPolytopes( ideal gls )
+pointSet ** convexHull::newtonPolytopesP( const ideal gls )
 {
   int i, j, k;
@@ -856 +841 @@
   return Q;
 }
+
+// mprSTICKYPROT:
+// ST_SPARSE_VADD: new vertex of convex hull added
+// ST_SPARSE_VREJ: point rejected (-> inside hull)
+ideal convexHull::newtonPolytopesI( const ideal gls )
+{
+  int i, j;
+  int m;  // Anzahl der Exponentvektoren im i-ten Polynom (gls->m)[i] des Ideals gls
+  int idelem= IDELEMS(gls);
+  ideal id;
+  poly p,pid,pd;
+  Exponent_t * vert;
+
+  n= pVariables;
+  vert= (Exponent_t *)Alloc( (idelem+1) * sizeof(Exponent_t) );
+  id= idInit( idelem, 1 );
+
+  for( i= 0; i < idelem; i++ )
+  {
+    m = pLength( (gls->m)[i] );
+
+    p= (gls->m)[i];
+    for( j= 1; j <= m; j++) {  // für jeden Exponentvektor
+      if( !inHull( (gls->m)[i], p, m, j ) )
+      {
+        if ( (id->m)[i] == NULL ) {
+          (id->m)[i]= pCopy(p);
+          pNext((id->m)[i])= NULL;
+          pid=(id->m)[i];
+        } else {
+          pNext(pid)= pCopy(p);
+          pIter(pid);
+          pNext(pid)= NULL;
+        }
+        mprSTICKYPROT(ST_SPARSE_VADD);
+      }
+      else
+      {
+        mprSTICKYPROT(ST_SPARSE_VREJ);
+      }
+      pIter( p );
+    } // j
+    mprSTICKYPROT("\n");
+  } // i
+
+  Free( (ADDRESS) vert, (idelem+1) * sizeof(Exponent_t) );
+
+#ifdef mprDEBUG_PROT
+  PrintLn();
+  for( i= 0; i < idelem; i++ )
+  {
+  }
+#endif
+
+  return id;
+}
 //<-
 
@@ -880 +921 @@
 {
   int i, ii, j, k, col, r;
-  int icase, constr, numverts, cols;
+  int numverts, cols;
 
   numverts = 0;
@@ -924 +965 @@
     Werror("mayanPyramidAlg::vDistance:"
            "setting up matrix for udist: col %d != cols %d",col,cols);
-  constr = n+dim+1;
+
+  pLP->m = n+dim+1;
+  pLP->m3= pLP->m;
+  pLP->n=cols-1;
 
 #ifdef mprDEBUG_ALL
-  Print("vDistance LP, known koords dim=%d, constr %d, cols %d, acoords= ",dim,constr,cols);
+  Print("vDistance LP, known koords dim=%d, constr %d, cols %d, acoords= ",dim,pLP->m,cols);
   for( i= 0; i < dim; i++ ) Print(" %d",acoords[i]);
   PrintLn();
-  print_mat( pLP->LiPM, constr+1, cols);
-#endif
-
-#if 1
-  if ( constr + 1 > pLP->LiPM_rows )
-    WerrorS("mayanPyramidAlg::vDistance: #rows > #pLP->LiPM_rows!");
-  if ( cols + 1 > pLP->LiPM_cols )
-    WerrorS("mayanPyramidAlg::vDistance: #cols > #pLP->LiPM_cols!");
-#endif
-
-  // 3rd arg is #columns excluding column of constants
-  //            M       N       m1  m2  m3,   M= m1+m2+m3
-  simplx( pLP->LiPM, constr, cols-1,  0,  0, constr, &icase, pLP->izrov, pLP->iposv);
+  print_mat( pLP->LiPM, pLP->m+1, cols);
+#endif
+  
+  pLP->compute();
 
 #ifdef mprDEBUG_ALL
   Print("LP returns matrix\n");
-  print_bmat( pLP->LiPM, constr+1, cols+1-constr, cols, pLP->iposv);
-#endif
-
-  if( icase != 0 ) {  // check for errors
+  print_bmat( pLP->LiPM, pLP->m+1, cols+1-pLP->m, cols, pLP->iposv);
+#endif
+
+  if( pLP->icase != 0 ) {  // check for errors
     WerrorS("mayanPyramidAlg::vDistance:");
-    if( icase == 1 )
-      WerrorS(" vDistance: Unbounded v-distance: probably 1st v-coor=0");
-    if( icase == -1 )
-      WerrorS(" vDistance: Infeasible v-distance");
+    if( pLP->icase == 1 )
+      WerrorS(" Unbounded v-distance: probably 1st v-coor=0");
+    else if( pLP->icase == -1 )
+      WerrorS(" Infeasible v-distance");
+    else
+      WerrorS(" Unknown error");
     return -1.0;
   }
@@ -964 +1001 @@
 {
   int i, j, k, cols, cons;
-  int icase, la_cons_row;
+  int la_cons_row;
 
   cons = n+dim+2;
@@ -1012 +1049 @@
 #endif
 
-#if 1
-  if ( cons + 1 > pLP->LiPM_rows )
-    WerrorS(" mn_mx_MinkowskiSum: #rows > #pLP->LiPM_rows!");
-  if ( cols + 1 > pLP->LiPM_cols )
-    WerrorS(" mn_mx_MinkowskiSum: #cols > #pLP->LiPM_cols!");
-#endif
-
   // simplx finds MIN for obj.fnc, puts it in [1,1]
-  simplx(pLP->LiPM, cons, cols-1, 0, 0, cons, &icase, pLP->izrov, pLP->iposv);
-
-  if ( icase != 0 ) { // check for errors
-    if( icase < 0)
+  pLP->m= cons;
+  pLP->n= cols-1;
+  pLP->m3= cons;
+
+  pLP->compute();
+
+  if ( pLP->icase != 0 ) { // check for errors
+    if( pLP->icase < 0)
       WerrorS(" mn_mx_MinkowskiSum: LinearProgram: minR: infeasible");
-    else if( icase > 0)
+    else if( pLP->icase > 0)
       WerrorS(" mn_mx_MinkowskiSum: LinearProgram: minR: unbounded");
   }
@@ -1071 +1105 @@
 #endif
 
-#if 1
-  if ( cons + 1 > pLP->LiPM_rows )
-    WerrorS(" mn_mx_MinkowskiSum: #rows > #pLP->LiPM_rows!");
-  if ( cols + 1 > pLP->LiPM_cols )
-    WerrorS(" mn_mx_MinkowskiSum: #cols > #pLP->LiPM_cols!");
-#endif
+  pLP->m= cons;
+  pLP->n= cols-1;
+  pLP->m3= cons;
 
   // simplx finds MAX for obj.fnc, puts it in [1,1]
-  simplx(pLP->LiPM, cons, cols-1, 0, 0, cons, &icase, pLP->izrov, pLP->iposv);
-
-  if ( icase != 0 )
-  {
-    if( icase < 0)
+  pLP->compute();
+
+  if ( pLP->icase != 0 )
+  {
+    if( pLP->icase < 0)
       WerrorS(" mn_mx_MinkowskiSum: LinearProgram: maxR: infeasible");
-    else if( icase > 0)
+    else if( pLP->icase > 0)
       WerrorS(" mn_mx_MinkowskiSum: LinearProgram: maxR: unbounded");
   }
@@ -1198 +1229 @@
 {
   int i, j, k,c ;
-  int ncols, size, icase, NumCons;
+  int size;
   bool found= true;
   mprfloat cd;
@@ -1205 +1236 @@
   setID *optSum;
 
-  ncols = 1;
-  NumCons = n + n + 1;   // number of constrains
+  LP->n = 1;
+  LP->m = n + n + 1;   // number of constrains
 
   // fill in LP matrix
@@ -1214 +1245 @@
     for ( k= 1; k <= size; k++ )
     {
-      ncols++;
+      LP->n++;
 
       // objective funtion, minimize
-      LP.LiPM[1][ncols] = - ( (mprfloat) (*pQ[i])[k]->point[pQ[i]->dim] / SCALEDOWN );
+      LP->LiPM[1][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[pQ[i]->dim] / SCALEDOWN );
 
       // lambdas sum up to 1
       for ( j = 0; j <= n; j++ )
         if ( i==j )
-          LP.LiPM[j+2][ncols] = -1.0;
+          LP->LiPM[j+2][LP->n] = -1.0;
         else
-          LP.LiPM[j+2][ncols] = 0.0;
+          LP->LiPM[j+2][LP->n] = 0.0;
 
       // the points
       for ( j = 1; j <= n; j++ )
       {
-        LP.LiPM[j+n+2][ncols] =  - ( (mprfloat) (*pQ[i])[k]->point[j] );
-      }
-    }
-  }
-
-  for ( j = 0; j <= n; j++ ) LP.LiPM[j+2][1] = 1.0;
+        LP->LiPM[j+n+2][LP->n] =  - ( (mprfloat) (*pQ[i])[k]->point[j] );
+      }
+    }
+  }
+
+  for ( j = 0; j <= n; j++ ) LP->LiPM[j+2][1] = 1.0;
   for ( j= 1; j <= n; j++ )
   {
-    LP.LiPM[j+n+2][1]= (mprfloat)(*E)[vert]->point[j] - shift[j];
-  }
-  ncols--;
-
-  LP.LiPM[1][1] = 0.0;
+    LP->LiPM[j+n+2][1]= (mprfloat)(*E)[vert]->point[j] - shift[j];
+  }
+  LP->n--;
+
+  LP->LiPM[1][1] = 0.0;
 
 #ifdef mprDEBUG_ALL
   PrintLn();
-  Print(" n= %d, NumCons=M= %d, ncols=N= %d\n",n,NumCons,ncols);
-  print_mat(LP.LiPM, NumCons+1, ncols+1);
-#endif
-
-#if 1
-  if ( NumCons + 1 > LP.LiPM_rows )
-    WerrorS("resMatrixSparse::RC: #rows > #LP.LiPM_rows!");
-  if ( ncols + 1 > LP.LiPM_cols )
-    WerrorS("resMatrixSparse::RC: #cols > #LP.LiPM_cols!");
-#endif
-
-  //           M        N      m1  m2  m3,   M= m1+m2+m3
-  simplx(LP.LiPM, NumCons, ncols,  0,  0, NumCons, &icase, LP.izrov, LP.iposv);
-
-  if (icase < 0)
+  Print(" n= %d, LP->m=M= %d, LP->n=N= %d\n",n,LP->m,LP->n);
+  print_mat(LP->LiPM, LP->m+1, LP->n+1);
+#endif
+
+  LP->m3= LP->m;
+
+  LP->compute();
+
+  if ( LP->icase < 0 )
   {
     // infeasibility: the point does not lie in a cell -> remove it
@@ -1266 +1291 @@
 
   // store result
-  (*E)[vert]->point[E->dim]= (int)(-LP.LiPM[1][1] * SCALEDOWN);
+  (*E)[vert]->point[E->dim]= (int)(-LP->LiPM[1][1] * SCALEDOWN);
 
 #ifdef mprDEBUG_ALL
-  Print(" simplx returned %d, Objective value = %f\n", icase, LP.LiPM[1][1]);
-  //print_bmat(LP.LiPM, NumCons + 1, ncols+1-NumCons, ncols+1, LP.iposv); // ( rows= M+1, cols= N+1-m3 )
-  //print_mat(LP.LiPM, NumCons+1, ncols);
+  Print(" simplx returned %d, Objective value = %f\n", LP->icase, LP->LiPM[1][1]);
+  //print_bmat(LP->LiPM, NumCons + 1, LP->n+1-NumCons, LP->n+1, LP->iposv); // ( rows= M+1, cols= N+1-m3 )
+  //print_mat(LP->LiPM, NumCons+1, LP->n);
 #endif
 
@@ -1279 +1304 @@
   {
     found=false;
-    for ( i= 1; i < NumCons; i++ )
-    {
-      if ( LP.iposv[i] > LP.iposv[i+1] )
-      {
-
-        c= LP.iposv[i];
-        LP.iposv[i]=LP.iposv[i+1];
-        LP.iposv[i+1]=c;
-
-        cd=LP.LiPM[i+1][1];
-        LP.LiPM[i+1][1]=LP.LiPM[i+2][1];
-        LP.LiPM[i+2][1]=cd;
+    for ( i= 1; i < LP->m; i++ )
+    {
+      if ( LP->iposv[i] > LP->iposv[i+1] )
+      {
+
+        c= LP->iposv[i];
+        LP->iposv[i]=LP->iposv[i+1];
+        LP->iposv[i+1]=c;
+
+        cd=LP->LiPM[i+1][1];
+        LP->LiPM[i+1][1]=LP->LiPM[i+2][1];
+        LP->LiPM[i+2][1]=cd;
 
         found= true;
@@ -1299 +1324 @@
 
 #ifdef mprDEBUG_ALL
-  print_bmat(LP.LiPM, NumCons + 1, ncols+1-NumCons, ncols+1, LP.iposv);
+  print_bmat(LP->LiPM, LP->m + 1, LP->n+1-LP->m, LP->n+1, LP->iposv);
   Print(" now split into sets\n");
 #endif
@@ -1308 +1333 @@
   // remap results of LP to sets Qi
   c=0;
-  optSum= (setID*)Alloc( (NumCons) * sizeof(struct setID) );
-  for ( i= 0; i < NumCons; i++ )
-  {
-    //Print("% .15f\n",LP.LiPM[i+2][1]);
-    if ( LP.LiPM[i+2][1] > 1e-12 )
-    {
-      if ( !remapXiToPoint( LP.iposv[i+1], pQ, &(optSum[c].set), &(optSum[c].pnt) ) )
-      {
-        Werror(" resMatrixSparse::RC: Found bad solution in LP: %d!",LP.iposv[i+1]);
+  optSum= (setID*)Alloc( (LP->m) * sizeof(struct setID) );
+  for ( i= 0; i < LP->m; i++ )
+  {
+    //Print("% .15f\n",LP->LiPM[i+2][1]);
+    if ( LP->LiPM[i+2][1] > 1e-12 )
+    {
+      if ( !remapXiToPoint( LP->iposv[i+1], pQ, &(optSum[c].set), &(optSum[c].pnt) ) )
+      {
+        Werror(" resMatrixSparse::RC: Found bad solution in LP: %d!",LP->iposv[i+1]);
         WerrorS(" resMatrixSparse::RC: remapXiToPoint faild!");
         return -1;
@@ -1348 +1373 @@
   if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
 
-#ifdef mprDEBUG_ALL
+#ifdef mprDEBUG_PROT
   Print("\n Point E[%d] was <",vert);print_exp((*E)[vert],E->dim-1);Print(">, bucket={");
   for ( j= 0; j < E->dim; j++ )
@@ -1355 +1380 @@
   }
   Print(" }\n optimal Sum: Qi ");
-  for ( j= 0; j < NumCons; j++ )
+  for ( j= 0; j < LP->m; j++ )
   {
     Print(" [ %d, %d ]",optSum[j].set,optSum[j].pnt);
@@ -1363 +1388 @@
 
   // clean up
-  Free( (ADDRESS) optSum, (NumCons) * sizeof(struct setID) );
+  Free( (ADDRESS) optSum, (LP->m) * sizeof(struct setID) );
 
   mprSTICKYPROT(ST_SPARSE_RC);
 
-  return (int)(-LP.LiPM[1][1] * SCALEDOWN);
+  return (int)(-LP->LiPM[1][1] * SCALEDOWN);
 }
 
@@ -1465 +1490 @@
   int i,j;
   i= 1;
-  /*
-  shift[1]= (mprfloat) (RVMULT*(siRand()%MAXRVVAL)/(mprfloat)MAXRVVAL);
-  i++;
-  while ( i <= dim )
-  {
-    shift[i]=shift[1];
-    i++;
-  }
-  */
+  time_t *tp = NULL;
+
   while ( i <= dim )
   {
@@ -1487 +1505 @@
     }
   }
-
 }
 
@@ -1567 +1584 @@
   idelem= IDELEMS(gls);  // should be n+1
 
-  // prepare matrix LP.LiPM for Linear Programming
+  // prepare matrix LP->LiPM for Linear Programming
   totverts = 0;
   for( i=0; i < idelem; i++) totverts += pLength( (gls->m)[i] );
 
-  LP.LiPM_rows= idelem+totverts*2+5;  // very approximal
-  LP.LiPM_cols= totverts+5;
-
-#ifdef mprDEBUG_ALL
-  Print("LP.LiPM size: %d, %d\n",LP.LiPM_rows,LP.LiPM_cols);
-#endif
-
-  // need AllocAligned since we allocate mem for type double
-  LP.LiPM = (mprfloat **)Alloc( LP.LiPM_rows * sizeof(mprfloat *) );  // LP matrix
-  for( i= 0; i < LP.LiPM_rows; i++ )
-  {
-    // Mem must be allocated aligned, also for type double!
-    LP.LiPM[i] = (mprfloat *)AllocAligned0( LP.LiPM_cols * sizeof(mprfloat) );
-    // LP.LiPM[i] = (mprfloat *)Alloc0( LP.LiPM_cols * sizeof(mprfloat) );
-  }
-
-  LP.iposv = (int *)Alloc0( (idelem * MAXPOINTS) * sizeof(int) );
-  LP.izrov = (int *)Alloc0( (idelem * MAXPOINTS) * sizeof(int) );
+  LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
 
   // get shift vector
@@ -1604 +1604 @@
 
   // evaluate convex hull for supports of gls
-  convexHull chnp( &LP );
-  Qi= chnp.newtonPolytopes( gls );
+  convexHull chnp( LP );
+  Qi= chnp.newtonPolytopesP( gls );
 
 #ifdef mprMINKSUM
@@ -1611 +1611 @@
 #else
   // get inner points
-  mayanPyramidAlg mpa( &LP );
+  mayanPyramidAlg mpa( LP );
   E= mpa.getInnerPoints( Qi, shift );
 #endif
 
-#ifdef mprDEBUG_ALL
+#ifdef mprDEBUG_PROT
 #ifdef mprMINKSUM
   Print("(MinkSum)");
@@ -1675 +1675 @@
   E->sort();
 
-#ifdef mprDEBUG_ALL
+#ifdef mprDEBUG_PROT
   Print(" points with a[ij] (%d):\n",E->num);
   for ( pnt= 1; pnt <= E->num; pnt++ )
@@ -1706 +1706 @@
   delete E;
 
-  for( i= 0; i < LP.LiPM_rows; i++ )
-  {
-    FreeAligned( (ADDRESS) LP.LiPM[i], LP.LiPM_cols * sizeof(mprfloat) );
-    //    Free( (ADDRESS) LP.LiPM[i], LP.LiPM_cols * sizeof(mprfloat) );
-  }
-  Free( (ADDRESS) LP.LiPM, LP.LiPM_rows * sizeof(mprfloat *) );
-
-  Free( (ADDRESS) LP.iposv, (idelem * MAXPOINTS) * sizeof(int) );
-  Free( (ADDRESS) LP.izrov, (idelem * MAXPOINTS) * sizeof(int) );
+  delete LP;
 }
 
@@ -2534 +2526 @@
 
   // evaluate determinant of matrix m using factory singclap_det
-  #ifdef HAVE_FACTORY
   poly res= singclap_det( m );
-  #else
-  poly res= mpDetBareiss( m );
-  #endif
 
   // avoid errors for det==0
@@ -2598 +2586 @@
   }
 
-  #ifdef HAVE_FACTORY
   poly res= singclap_det( mat );
-  #else
-  poly res= mpDetBareiss( mat );
-  #endif
 
   number numres;
@@ -2623 +2607 @@
 
 #define MAXEVPOINT 1000000 // 0x7fffffff
+//#define MPR_MASI
 
 //-> unsigned long over(unsigned long n,unsigned long d)
@@ -2938 +2923 @@
         else if ( i <= uvar + 2 )
         {
-          pevpoint[i]=nInit(1+(siRand() % MAXEVPOINT));
+          pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
           //pevpoint[i]=nInit(383);
         }
@@ -3058 +3043 @@
         if ( i <= uvar + 2 )
         {
-          pevpoint[i]=nInit(1+(siRand() % MAXEVPOINT));
+          pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
           //pevpoint[i]=nInit(383);
         } else pevpoint[i]=nInit(0);
@@ -3080 +3065 @@
     number *ncpoly= (number *)Alloc( (tdg+1) * sizeof(number) );
 
-#ifdef MPR_TIMING
-       BOOLEAN masi=true;
-#endif
+#ifdef MPR_MASI
+    BOOLEAN masi=true;
+#endif
+
     piter= pures;
     for ( i= tdg; i >= 0; i-- )
@@ -3091 +3077 @@
         ncpoly[i]= nCopy( pGetCoeff( piter ) );
         pIter( piter );
-#ifdef MPR_TIMING
-       masi=false;
+#ifdef MPR_MASI
+        masi=false;
 #endif
       }
@@ -3101 +3087 @@
       mprPROTNnl("", ncpoly[i] );
     }
-
-#ifdef MPR_TIMING
-     if ( masi ) Print("MASI MASI MASI\n");
+#ifdef MPR_MASI
+    if ( masi ) mprSTICKYPROT("MASI");
 #endif
 
@@ -3149 +3134 @@
 }
 //<-
+
+//-----------------------------------------------------------------------------
+
+//-> loNewtonPolytope(...)
+ideal loNewtonPolytope( const ideal id )
+{
+  simplex * LP;
+  int i;
+  int n,totverts,idelem;
+  ideal idr;
+
+  n= pVariables;
+  idelem= IDELEMS(id);  // should be n+1
+
+  totverts = 0;
+  for( i=0; i < idelem; i++) totverts += pLength( (id->m)[i] );
+
+  LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
+
+  // evaluate convex hull for supports of id
+  convexHull chnp( LP );
+  idr = chnp.newtonPolytopesI( id );
+
+  delete LP;
+
+  return idr;
+}
+//<-
+
 //%e
 
@@ -3164 +3178 @@
 // leave fold: C-c y
 //   foldmode: F10
-

Singular/mpr_base.h

@@ -4 +4 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: mpr_base.h,v 1.4 1999-11-15 17:20:30 obachman Exp $ */
+/* $Id: mpr_base.h,v 1.5 1999-12-02 23:03:50 wenk Exp $ */
 
 /*
@@ -94 +94 @@
 };
 //<-
+
+ideal loNewtonPolytope( const ideal id );
 //%e
 #endif MPR_BASE_H

Singular/mpr_global.h

@@ -4 +4 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: mpr_global.h,v 1.8 1999-11-15 17:20:30 obachman Exp $ */
+/* $Id: mpr_global.h,v 1.9 1999-12-02 23:03:50 wenk Exp $ */
 
 /*
@@ -12 +12 @@
 
 // to get detailed timigs, define MPR_TIMING
-// #define MPR_TIMING
+//#define MPR_TIMING
 
 // Set to double or long double. double is recomended.

Singular/mpr_inout.cc

@@ -3 +3 @@
 ****************************************/
 
-/* $Id: mpr_inout.cc,v 1.6 1999-11-15 17:20:30 obachman Exp $ */
+/* $Id: mpr_inout.cc,v 1.7 1999-12-02 23:03:51 wenk Exp $ */
 
 /*
@@ -25 +25 @@
 #include "numbers.h"
 #include "lists.h"
+#include "matpol.h"
 
 #include <math.h>
@@ -36 +37 @@
 #ifdef MPR_TIMING
 #define TIMING
+#endif
 #include "../factory/timing.h"
 TIMING_DEFINE_PRINT(mpr_overall)
@@ -44 +46 @@
 TIMING_DEFINE_PRINT(mpr_arrange)
 TIMING_DEFINE_PRINT(mpr_solver)
+
 #define TIMING_EPR(t,msg) TIMING_END_AND_PRINT(t,msg);TIMING_RESET(t);
-#else
-#define TIMING_START(x)
-#define TIMING_EPR(x,y)
-#endif
-
 
 enum mprState
@@ -181 +179 @@
   if ( v->Typ() != IDEAL_CMD )
     return TRUE;
-  else gls= (ideal)(args->Data());
+  else gls= (ideal)(v->Data());
   v= v->next;
 
@@ -211 +209 @@
   BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
 
-  //emptylist= (lists)AllocSizeOf( slists );
+  //emptylist= (lists)Alloc( sizeof(slists) );
   //emptylist->Init( 0 );
 
@@ -314 +312 @@
 
   //emptylist->Clean();
-  //  FreeSizeOf( (ADDRESS) emptylist, slists );
+  //  Free( (ADDRESS) emptylist, sizeof(slists) );
 
   TIMING_EPR(mpr_overall,"overall time\t\t")
@@ -356 +354 @@
   int howclean= (int)arg3->Data();
 
+  if ( !(rField_is_R() ||
+         rField_is_Q() ||
+         rField_is_long_R() ||
+         rField_is_long_C()) )
+  {
+    WerrorS("Ground field not implemented!");
+    return TRUE;
+  }
+  
   if ( !(rField_is_R()||rField_is_long_R()||rField_is_long_C()) ) 
   {
     setGMPFloatDigits( (unsigned long int)arg2->Data() );
+  }
+
+  if ( gls == NULL || pIsConstant( gls ) )
+  {
+    WerrorS("Input polynomial is constant!");
+    return TRUE;
   }
 
@@ -369 +382 @@
   lists rlist;
 
-  elist= (lists)AllocSizeOf( slists );
+  elist= (lists)Alloc( sizeof(slists) );
   elist->Init( 0 );
-
-  if ( !(rField_is_R() ||
-         rField_is_Q() ||
-         rField_is_long_R() ||
-         rField_is_long_C()) )
-         {
-    WerrorS("Ground field not implemented!");
-    return TRUE;
-  }
 
   if ( pVariables > 1 )
@@ -395 +399 @@
         if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
         {
-          Werror("The polynomial %s must be univariate!",arg1->Name());
+          WerrorS("The input polynomial must be univariate!");
           return TRUE;
         }
@@ -438 +442 @@
   int j;
 
-  rlist= (lists)AllocSizeOf( slists );
+  rlist= (lists)Alloc( sizeof(slists) );
   rlist->Init( elem );
 
@@ -461 +465 @@
 
   elist->Clean();
-  //FreeSizeOf( (ADDRESS) elist, slists );
+  //Free( (ADDRESS) elist, sizeof(slists) );
 
   for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
@@ -629 +633 @@
 //<-
 
+//-> BOOLEAN loNewtonP( leftv res, leftv arg1 )
+BOOLEAN loNewtonP( leftv res, leftv arg1 )
+{
+  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
+  return FALSE;
+}
+//<-
+
+//-> BOOLEAN loSimplex( leftv res, leftv args )
+BOOLEAN loSimplex( leftv res, leftv args )
+{
+  if ( !(rField_is_long_R()) )
+  {
+    WerrorS("Ground field not implemented!");
+    return TRUE;
+  }
+
+  simplex * LP;
+  matrix m;
+
+  leftv v= args;
+  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
+    return TRUE;
+  else 
+    m= (matrix)(v->Data());
+  
+  LP = new simplex(MATROWS(m),MATCOLS(m));
+  LP->mapFromMatrix(m);
+
+  v= v->next;
+  if ( v->Typ() != INT_CMD )    // 2: m = number of constraints
+    return TRUE;
+  else 
+    LP->m= (int)(v->Data());
+
+  v= v->next;
+  if ( v->Typ() != INT_CMD )    // 3: n = number of variables 
+    return TRUE;
+  else 
+    LP->n= (int)(v->Data());
+
+  v= v->next;
+  if ( v->Typ() != INT_CMD )    // 4: m1 = number of <= constraints
+    return TRUE;
+  else 
+    LP->m1= (int)(v->Data());
+
+  v= v->next;
+  if ( v->Typ() != INT_CMD )    // 5: m2 = number of >= constraints
+    return TRUE;
+  else 
+    LP->m2= (int)(v->Data());
+
+  v= v->next;
+  if ( v->Typ() != INT_CMD )    // 6: m3 = number of == constraints
+    return TRUE;
+  else 
+    LP->m3= (int)(v->Data());
+
+#ifdef mprDEBUG_PROT
+  Print("m (constraints) %d\n",LP->m);
+  Print("n (columns) %d\n",LP->n);
+  Print("m1 (<=) %d\n",LP->m1);
+  Print("m2 (>=) %d\n",LP->m2);
+  Print("m3 (==) %d\n",LP->m3);
+#endif
+
+  LP->compute();
+
+  lists lres= (lists)Alloc( sizeof(slists) );
+  lres->Init( 6 );
+
+  lres->m[0].rtyp= MATRIX_CMD; // output matrix
+  lres->m[0].data=(void*)LP->mapToMatrix(m);
+
+  lres->m[1].rtyp= INT_CMD;   // found a solution?
+  lres->m[1].data=(void*)LP->icase;
+
+  lres->m[2].rtyp= INTVEC_CMD;
+  lres->m[2].data=(void*)LP->posvToIV();
+ 
+  lres->m[3].rtyp= INTVEC_CMD;
+  lres->m[3].data=(void*)LP->zrovToIV();
+
+  lres->m[4].rtyp= INT_CMD;
+  lres->m[4].data=(void*)LP->m;
+
+  lres->m[5].rtyp= INT_CMD;
+  lres->m[5].data=(void*)LP->n;
+
+  res->data= (void*)lres;
+
+  return FALSE;
+}
+//<-
+
 //-----------------------------------------------------------------------------
 

Singular/mpr_inout.h

@@ -5 +5 @@
 ****************************************/
 
-/* $Id: mpr_inout.h,v 1.5 1999-11-15 17:20:31 obachman Exp $ */
+/* $Id: mpr_inout.h,v 1.6 1999-12-02 23:03:51 wenk Exp $ */
 
 /*
@@ -58 +58 @@
 BOOLEAN nuVanderSys( leftv res, leftv arg1, leftv arg2, leftv arg3 );
 
+/** compute Newton Polytopes of input polynomials
+ */
+BOOLEAN loNewtonP( leftv res, leftv arg1 );
+
+/** Implementation of the Simplex Algorithm.
+ * For args, see class simplex.
+ */
+BOOLEAN loSimplex( leftv res, leftv args );
+
 #endif
 

Singular/mpr_numeric.cc

@@ -3 +3 @@
 ****************************************/
 
-/* $Id: mpr_numeric.cc,v 1.6 1999-11-15 17:20:31 obachman Exp $ */
+/* $Id: mpr_numeric.cc,v 1.7 1999-12-02 23:03:51 wenk Exp $ */
 
 /*
@@ -24 +24 @@
 #include "intvec.h"
 #include "longalg.h"
+#include "matpol.h"
 #include "ring.h"
 //#include "longrat.h"
@@ -749 +750 @@
   int elem= roots[0]->getAnzElems();  // number of koordinates per root
 
-  lists listofroots= (lists)AllocSizeOf( slists ); // must be done this way!
+  lists listofroots= (lists)Alloc( sizeof(slists) ); // must be done this way!
 
   if ( found_roots )
@@ -757 +758 @@
     for (i=0; i < count; i++)
     {
-      lists onepoint= (lists)AllocSizeOf(slists); // must be done this way!
+      lists onepoint= (lists)Alloc(sizeof(slists)); // must be done this way!
       onepoint->Init(elem);
       for ( j= 0; j < elem; j++ )
@@ -791 +792 @@
 
 //-----------------------------------------------------------------------------
-//-------------- ludcmp/lubksb ------------------------------------------------
+//-------------- simplex ----- ------------------------------------------------
 //-----------------------------------------------------------------------------
 
-//#define error(a) a
-#define error(a)
-
-//-> simplex
+//  #ifdef mprDEBUG_PROT
+//  #define error(a) a
+//  #else
+//  #define error(a)
+//  #endif
+
+#define error(a) a
+
+#define MAXPOINTS      1000
+
+//-> simplex::*
 //
-void simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax );
-void simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 );
-void simp3( mprfloat **a, int i1, int k1, int ip, int kp );
-
-void simplx( mprfloat **a, int m, int n, int m1, int m2, int m3, int *icase, int izrov[], int iposv[] )
+simplex::simplex( int rows, int cols )
+   : LiPM_cols(cols), LiPM_rows(rows)
+{
+  int i;
+
+  LiPM_rows=LiPM_rows+3;
+  LiPM_cols=LiPM_cols+2;
+
+  LiPM = (mprfloat **)Alloc( LiPM_rows * sizeof(mprfloat *) );  // LP matrix
+  for( i= 0; i < LiPM_rows; i++ )
+  {
+    // Mem must be allocated aligned, also for type double!
+    LiPM[i] = (mprfloat *)AllocAligned0( LiPM_cols * sizeof(mprfloat) );
+  }
+
+  iposv = (int *)Alloc0( 2*LiPM_rows*sizeof(int) );
+  izrov = (int *)Alloc0( 2*LiPM_rows*sizeof(int) );
+
+  m=n=m1=m2=m3=icase=0;
+
+#ifdef mprDEBUG_ALL
+  Print("LiPM size: %d, %d\n",LiPM_rows,LiPM_cols);
+#endif
+}
+
+simplex::~simplex()
+{
+  // clean up
+  int i;
+  for( i= 0; i < LiPM_rows; i++ )
+  {
+    FreeAligned( (ADDRESS) LiPM[i], LiPM_cols * sizeof(mprfloat) );
+  }
+  Free( (ADDRESS) LiPM, LiPM_rows * sizeof(mprfloat *) );
+
+  Free( (ADDRESS) iposv, 2*LiPM_rows*sizeof(int) );
+  Free( (ADDRESS) izrov, 2*LiPM_rows*sizeof(int) );
+}
+
+BOOLEAN simplex::mapFromMatrix( matrix m )
+{
+  int i,j;
+//    if ( MATROWS( m ) > LiPM_rows ||  MATCOLS( m ) > LiPM_cols ) {
+//      WerrorS("");
+//      return FALSE;
+//    }
+
+  number coef;
+  for ( i= 1; i <= MATROWS( m ); i++ )
+  {
+     for ( j= 1; j <= MATCOLS( m ); j++ ) 
+     {
+        if ( MATELEM(m,i,j) != NULL ) 
+        {
+           coef= pGetCoeff( MATELEM(m,i,j) );
+           if ( coef != NULL && !nIsZero(coef) )
+              LiPM[i][j]= (double)(*(gmp_float*)coef);
+           //#ifdef mpr_DEBUG_PROT
+           //Print("%f ",LiPM[i][j]);
+           //#endif
+        }
+     }
+     //     PrintLn();
+  }
+
+  return TRUE;
+}
+
+matrix simplex::mapToMatrix( matrix m )
+{
+  int i,j;
+//    if ( MATROWS( m ) < LiPM_rows-3 ||  MATCOLS( m ) < LiPM_cols-2 ) {
+//      WerrorS("");
+//      return NULL;
+//    }
+
+//Print(" %d x %d\n",MATROWS( m ),MATCOLS( m ));
+
+  number coef;
+  gmp_float * bla;
+  for ( i= 1; i <= MATROWS( m ); i++ )
+  {
+    for ( j= 1; j <= MATCOLS( m ); j++ )
+    {
+       pDelete( &(MATELEM(m,i,j)) );
+       MATELEM(m,i,j)= NULL;
+//Print(" %3.0f ",LiPM[i][j]);
+       if ( LiPM[i][j] != 0.0 ) 
+       {
+          bla= new gmp_float(LiPM[i][j]);
+          coef= (number)bla;
+          MATELEM(m,i,j)= pOne();
+          pSetCoeff( MATELEM(m,i,j), coef );
+       }
+    }
+//PrintLn();
+  }
+  
+  return m;
+}
+
+intvec * simplex::posvToIV()
+{
+   int i;
+   intvec * iv = new intvec( m );
+   for ( i= 1; i <= m; i++ )
+   {
+      IMATELEM(*iv,i,1)= iposv[i];
+   }
+   return iv;
+}
+
+intvec * simplex::zrovToIV()
+{
+   int i;
+   intvec * iv = new intvec( n );
+   for ( i= 1; i <= n; i++ )
+   {
+      IMATELEM(*iv,i,1)= izrov[i];
+   }
+   return iv; 
+}
+
+void simplex::compute()
 {
   int i,ip,ir,is,k,kh,kp,m12,nl1,nl2;
@@ -809 +936 @@
   mprfloat q1, bmax;
 
-  if ( m != (m1+m2+m3))
+  if ( m != (m1+m2+m3) )
   {
     // error: bad input
-    error(WerrorS(" bad input constraint counts in simplex ");)
-    *icase=-2;
+    error(WerrorS("simplex::compute: Bad input constraint counts!");)
+    icase=-2;
     return;
   }
@@ -826 +953 @@
   for ( i=1; i<=m; i++ )
   {
-    if ( a[i+1][1] < 0.0 )
+    if ( LiPM[i+1][1] < 0.0 )
     {
       // error: bad input
-      error(WerrorS(" bad input tableau in simplex ");)
-      *icase=-2;
+      error(WerrorS("simplex::compute: Bad input tableau!");)
+      error(Werror("simplex::compute: in input Matrix row %d, column 1, value %f",i+1,LiPM[i+1][1]);)
+      icase=-2;
       // free mem l1,l2,l3;
       Free( (ADDRESS) l3, (m+1) * sizeof(int) );
@@ -848 +976 @@
     {
       q1=0.0;
-      for ( i=m1+1; i <= m; i++ ) q1+= a[i+1][k];
-      a[m+2][k]= -q1;
+      for ( i=m1+1; i <= m; i++ ) q1+= LiPM[i+1][k];
+      LiPM[m+2][k]= -q1;
     }
 
     do
     {
-      simp1(a,m+1,l1,nl1,0,&kp,&bmax);
-      if ( bmax <= SIMPLEX_EPS && a[m+2][1] < -SIMPLEX_EPS )
-      {
-        *icase= -1; // no solution found
+      simp1(LiPM,m+1,l1,nl1,0,&kp,&bmax);
+      if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] < -SIMPLEX_EPS )
+      {
+        icase= -1; // no solution found
         // free mem l1,l2,l3;
         Free( (ADDRESS) l3, (m+1) * sizeof(int) );
@@ -864 +992 @@
         return;
       }
-      else if ( bmax <= SIMPLEX_EPS && a[m+2][1] <= SIMPLEX_EPS )
+      else if ( bmax <= SIMPLEX_EPS && LiPM[m+2][1] <= SIMPLEX_EPS )
       {
         m12= m1+m2+1;
@@ -873 +1001 @@
             if ( iposv[ip] == (ip+n) )
             {
-              simp1(a,ip,l1,nl1,1,&kp,&bmax);
+              simp1(LiPM,ip,l1,nl1,1,&kp,&bmax);
               if ( fabs(bmax) >= SIMPLEX_EPS)
                 goto one;
@@ -885 +1013 @@
             if ( l3[i-m1] == 1 )
               for ( k=1; k <= n+1; k++ )
-                a[i+1][k] = -a[i+1][k];
+                LiPM[i+1][k] = -(LiPM[i+1][k]);
         break;
       }
@@ -891 +1019 @@
       //print_bmat( a, m+2, n+3);
       //#endif
-      simp2(a,n,l2,nl2,&ip,kp,&q1);
+      simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
       if ( ip == 0 )
       {
-        *icase = -1; // no solution found
+        icase = -1; // no solution found
         // free mem l1,l2,l3;
         Free( (ADDRESS) l3, (m+1) * sizeof(int) );
@@ -901 +1029 @@
         return;
       }
-    one: simp3(a,m+1,n,ip,kp);
+    one: simp3(LiPM,m+1,n,ip,kp);
       // #if DEBUG
       // print_bmat(a,m+2,n+3);
@@ -907 +1035 @@
       if ( iposv[ip] >= (n+m1+m2+1))
       {
-        for ( k=1; k<= nl1; k++ )
+        for ( k= 1; k <= nl1; k++ )
           if ( l1[k] == kp ) break;
         --nl1;
         for ( is=k; is <= nl1; is++ ) l1[is]= l1[is+1];
-        ++a[m+2][kp+1];
-        for ( i= 1; i <= m+2; i++ ) a[i][kp+1] = -a[i][kp+1];
+        ++(LiPM[m+2][kp+1]);
+        for ( i= 1; i <= m+2; i++ ) LiPM[i][kp+1] = -(LiPM[i][kp+1]);
       }
       else
@@ -922 +1050 @@
           {
             l3[kh]= 0;
-            ++a[m+2][kp+1];
+            ++(LiPM[m+2][kp+1]);
             for ( i=1; i<= m+2; i++ )
-              a[i][kp+1] = -a[i][kp+1];
+              LiPM[i][kp+1] = -(LiPM[i][kp+1]);
           }
         }
@@ -939 +1067 @@
     // print_bmat( a, m+1, n+5);
     // #endif
-    simp1(a,0,l1,nl1,0,&kp,&bmax);
+    simp1(LiPM,0,l1,nl1,0,&kp,&bmax);
     if (bmax <= /*SIMPLEX_EPS*/0.0)
     {
-      *icase=0; // finite solution found
+      icase=0; // finite solution found
       // free mem l1,l2,l3
       Free( (ADDRESS) l3, (m+1) * sizeof(int) );
@@ -949 +1077 @@
       return;
     }
-    simp2(a,n,l2,nl2,&ip,kp,&q1);
+    simp2(LiPM,n,l2,nl2,&ip,kp,&q1);
     if (ip == 0)
     {
@@ -956 +1084 @@
       //       print_bmat( a, m+1, n+1);
       // #endif
-      *icase=1;                /* unbounded */
+      icase=1;                /* unbounded */
       // free mem
       Free( (ADDRESS) l3, (m+1) * sizeof(int) );
@@ -963 +1091 @@
       return;
     }
-    simp3(a,m,n,ip,kp);
+    simp3(LiPM,m,n,ip,kp);
     is= izrov[kp];
     izrov[kp]= iposv[ip];
@@ -970 +1098 @@
 }
 
-void simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax )
+void simplex::simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax )
 {
   int k;
@@ -1005 +1133 @@
 }
 
-void simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 )
+void simplex::simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 )
 {
   int k,ii,i;
@@ -1044 +1172 @@
 }
 
-void simp3( mprfloat **a, int i1, int k1, int ip, int kp )
+void simplex::simp3( mprfloat **a, int i1, int k1, int ip, int kp )
 {
   int kk,ii;

Singular/mpr_numeric.h

@@ -5 +5 @@
 ****************************************/
 
-/* $Id: mpr_numeric.h,v 1.4 1999-11-15 17:20:31 obachman Exp $ */
+/* $Id: mpr_numeric.h,v 1.5 1999-12-02 23:03:52 wenk Exp $ */
 
 /*
@@ -162 +162 @@
 #define SIMPLEX_EPS 1.0e-12
 
-void simplx( mprfloat **a, int m, int n, int m1, int m2, int m3, int *icase, int izrov[], int iposv[] );
+/** Linear Programming / Linear Optimization using Simplex - Algorithm
+ *
+ * On output, the tableau LiPM is indexed by two arrays of integers.
+ * ipsov[j] contains, for j=1..m, the number i whose original variable
+ * is now represented by row j+1 of LiPM. (left-handed vars in solution)
+ * (first row is the one with the objective function)
+ * izrov[j] contains, for j=1..n, the number i whose original variable
+ * x_i is now a right-handed variable, rep. by column j+1 of LiPM.
+ * These vars are all zero in the solution. The meaning of n<i<n+m1+m2
+ * is the same as above.
+ */
+class simplex
+{
+public:
+  
+  int m;         // number of constraints, make sure m == m1 + m2 + m3 !!
+  int n;         // # of independent variables
+  int m1,m2,m3;  // constraints <=, >= and ==
+  int icase;     // == 0: finite solution found; 
+                 // == +1 objective funtion unbound; == -1: no solution
+  int *izrov,*iposv;    
+
+  mprfloat **LiPM; // the matrix (of size [m+2, n+1])
+
+  /** #rows should be >= m+2, #cols >= n+1 
+   */
+  simplex( int rows, int cols );
+  ~simplex();
+
+  BOOLEAN mapFromMatrix( matrix m );
+  matrix mapToMatrix( matrix m );
+  intvec * posvToIV();
+  intvec * zrovToIV();
+
+  void compute();
+
+private:
+  simplex( const simplex & );
+  void simp1( mprfloat **a, int mm, int ll[], int nll, int iabf, int *kp, mprfloat *bmax );
+  void simp2( mprfloat **a, int n, int l2[], int nl2, int *ip, int kp, mprfloat *q1 );
+  void simp3( mprfloat **a, int i1, int k1, int ip, int kp );
+
+  int LiPM_cols,LiPM_rows;
+};
+
 //<-
 

Singular/tok.h

@@ -7 +7 @@
 * ABSTRACT: tokens, types for interpreter; general macros
 */
-/* $Id: tok.h,v 1.33 1999-11-29 14:46:55 Singular Exp $ */
+/* $Id: tok.h,v 1.34 1999-12-02 23:03:52 wenk Exp $ */
 
 #ifndef MYYSTYPE
@@ -89 +89 @@
   NAMEOF_CMD,
   NAMES_CMD,
+  NEWTONPOLY_CMD,
   NPARS_CMD,
   NVARS_CMD,
@@ -106 +107 @@
   RESULTANT_CMD,
   ROWS_CMD,
+  SIMPLEX_CMD,
   SQR_FREE_DEC_CMD,
   STATUS_CMD,