Changeset e858e7 – Singular

Singular/clapsing.cc

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 // $Id: clapsing.cc,v 1.52 1999-06-15 08:32:31 Singular Exp $
+// $Id: clapsing.cc,v 1.53 1999-06-28 16:06:22 Singular Exp $
 /*
 * ABSTRACT: interface between Singular and factory
 …
+    {
       (*v)=new intvec(1);
       (*v)[1]=1;
+      (**v)[0]=1;
+    }
     return res;
 …
   ideal f=singclap_factorize((poly)(u->Data()), &v, 0);
   if (f==NULL) return TRUE;
+  #ifdef MDEBUG
+  v->ivTEST();
+  #endif
   lists l=(lists)Alloc(sizeof(slists));
   l->Init(2);

Singular/intvec.cc

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR      *
 *****************************************/
 /* $Id: intvec.cc,v 1.12 1999-04-17 12:30:17 Singular Exp $ */
+/* $Id: intvec.cc,v 1.13 1999-06-28 16:06:23 Singular Exp $ */
 /*
 * ABSTRACT: class intvec: lists/vectors of integers
 …
   if ((col == 1)&&(mat!=INTMAT_CMD))
+  {
     int i;
     for (i=0; i<row-1; i++)
+    int i=0;
+    for (; i<row-1; i++)
+    {
       StringAppend("%d,",v[i]);

Singular/intvec.h

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: intvec.h,v 1.7 1999-04-16 07:53:34 obachman Exp $ */
+/* $Id: intvec.h,v 1.8 1999-06-28 16:06:23 Singular Exp $ */
 /*
 * ABSTRACT: class intvec: lists/vectors of integers
 …
     int range(int i, int j)
        { return ((i<row) && (i>=0) && (j<col) && (j>=0)); }
+    int& operator[](int i) {
+    int& operator[](int i)
+    {
+       #ifndef NDEBUG
        if((i<0)||(i>=row*col))
+       {
           Werror("wrong intvec index:%d\n",i);
+       }
+       return v[i];}
+       #endif
+       return v[i];
+    }
 #define IMATELEM(M,I,J) (M)[(I-1)*(M).cols()+J-1]
     void operator+=(int intop);
 …
+         }
+       }
+    void ivTEST()
+       {
+#ifdef MDEBUG
+         mmTestL(this);
+         mmTest((ADDRESS)v,sizeof(int)*row*col);
+#endif
+       }
 };
 intvec * ivCopy(intvec * o);

Singular/mpr_base.cc

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: mpr_base.cc,v 1.1 1999-06-28 12:48:11 wenk Exp $ */
 /*
+/* $Id: mpr_base.cc,v 1.2 1999-06-28 16:06:24 Singular Exp $ */
+/*
  * ABSTRACT - multipolynomial resultants - resultant matrices
  *            ( sparse, dense, u-resultant solver )
 …
 #define MINVDIST       0.0
 #define RVMULT         0.00005 // 0.000000005
 #define MAXVARS        100
+#define MAXVARS        100
 //<-
 …
 /* Linear Program stuff */
+struct linProg {          // Linear Program stuff
+  mprfloat **LiPM;
+  int *izrov, *iposv;
+  int LiPM_cols,LiPM_rows;
+struct linProg
+{
+  mprfloat **LiPM;
+  int *izrov, *iposv;
+  int LiPM_cols,LiPM_rows;
 };
 …
   const poly getUDet( const number* evpoint );
 private:
   resMatrixSparse( const resMatrixSparse & );
 …
   /* Remaps a result of LP to the according point set Qi.
    * Returns false iff remaping was not possible, otherwise true.
+   * Returns false iff remaping was not possible, otherwise true.
    */
   bool remapXiToPoint( const int indx, pointSet **pQ, int *set, int *vtx );
 …
 private:
   ideal gls;
   int n, idelem;     // number of variables, polynoms
   int numSet0;       // number of elements in S0
   int msize;         // size of matrix
   intvec *uRPos;
   ideal rmat;        // sparse matrix representation
   linProg LP;
 };
 …
 typedef unsigned int Coord_t;
+struct setID {
+struct setID
+{
   int set;
   int pnt;
 };
+struct onePoint {
+struct onePoint
+{
   Coord_t * point;             // point[0] is unused, maxial dimension is MAXVARS+1
   setID rc;                    // filled in by Row Content Function
 …
 typedef struct onePoint * onePointP;
+/* sparse matrix entry */
+struct _entry {
+/* sparse matrix entry */
+struct _entry
+{
   number num;
   int col;
 …
 //-> class pointSet
+class pointSet {
+class pointSet
+{
 private:
   onePointP *points;     // set of onePoint's, index [1..num], supports of monoms
 …
   /* Returns index of supp(LT(p)) in pointSet. */
   int getExpPos( const poly p );
   /** sort lex
    */
 …
 //-> class convexHull
 /* Compute convex hull of given exponent set */
+class convexHull {
+class convexHull
+{
 public:
   convexHull( linProgP _pLP ) : pLP(_pLP) {}
 …
 private:
   /** Returns true iff the support of poly pointPoly is inside the
+  /** Returns true iff the support of poly pointPoly is inside the
    * convex hull of all points given by the  support of poly p.
    */
 …
 //-> class mayanPyramidAlg
 /* Compute all lattice points in a given convex hull */
+class mayanPyramidAlg {
+class mayanPyramidAlg
+{
 public:
   mayanPyramidAlg( linProgP _pLP ) : n(pVariables), pLP(_pLP) {}
 …
 private:
   /** Recursive Mayan Pyramid algorithm for directly computing MinkowskiSum
    * lattice points for (n+1)-fold MinkowskiSum of given point sets Qi[].
+  /** Recursive Mayan Pyramid algorithm for directly computing MinkowskiSum
+   * lattice points for (n+1)-fold MinkowskiSum of given point sets Qi[].
    * Recursively for range of dim: dim in [0..n); acoords[0..var) fixed.
    * Stores only MinkowskiSum points of udist > 0: done by storeMinkowskiSumPoints.
 …
   /**  Compute v-distance via Linear Programing
    * Linear Program finds the v-distance of the point in accords[].
    * The v-distance is the distance along the direction v to boundary of
+   * The v-distance is the distance along the direction v to boundary of
    * Minkowski Sum of Qi (here vector v is represented by shift[]).
    * Returns the v-distance or -1.0 if an error occured.
 …
   /** LP for finding min/max coord in MinkowskiSum, given previous coors.
    * Assume MinkowskiSum in non-negative quadrants
    * coor in [0,n); fixed coords in acoords[0..coor)
+   * coor in [0,n); fixed coords in acoords[0..coor)
    */
   void mn_mx_MinkowskiSum( int dim, Coord_t *minR, Coord_t *maxR );
 …
    */
   bool storeMinkowskiSumPoint();
 private:
   pointSet **Qi;
 …
   int i, j;
+  for (i = 1; i <= maxrow; i++) {
+  for (i = 1; i <= maxrow; i++)
+  {
     Print("[");
     for (j = 1; j <= maxcol; j++) Print("% 7.2f, ", a[i][j]);
 …
   printf("Output matrix from LinProg");
+  for (i = 1; i <= nrows; i++) {
+  for (i = 1; i <= nrows; i++)
+  {
     printf("\n[ ");
     if (i == 1) printf("  ");
 …
+{
   int i;
+  for ( i= 1; i <= n; i++ ) {
+  for ( i= 1; i <= n; i++ )
+  {
     Print(" %d",vert->point[i] );
 #ifdef LONG_OUTPUT
 …
   int val;
   Print(" matrix m[%d][%d]=(\n",MATROWS( omat ),MATCOLS( omat ));
+  for ( i= 1; i <= MATROWS( omat ); i++ ) {
+    for ( j= 1; j <= MATCOLS( omat ); j++ ) {
+      if ( MATELEM( omat, i, j) && pGetCoeff( MATELEM( omat, i, j) ) ) {
+        val= nInt(pGetCoeff( MATELEM( omat, i, j) ));
+        if ( i==MATROWS(omat) && j==MATCOLS(omat) ) {
+          if ( !val ) Print(" 0 ");
+          else Print("%d ",val);
+        } else {
+          if ( !val ) Print(" 0, ");
+          else Print("%d, ",val);
+        }
+      } else {
+        if ( i==MATROWS(omat) && j==MATCOLS(omat) ) {
+          Print("  0");
+        } else {
+          Print("  0, ");
+        }
+  for ( i= 1; i <= MATROWS( omat ); i++ )
+  {
+    for ( j= 1; j <= MATCOLS( omat ); j++ )
+    {
+      if ( MATELEM( omat, i, j) && pGetCoeff( MATELEM( omat, i, j) ) )
+      {
+        val= nInt(pGetCoeff( MATELEM( omat, i, j) ));
+        if ( i==MATROWS(omat) && j==MATCOLS(omat) )
+        {
+          if ( !val ) Print(" 0 ");
+          else Print("%d ",val);
+        }
+        else
+        {
+          if ( !val ) Print(" 0, ");
+          else Print("%d, ",val);
+        }
+      }
+      else
+      {
+        if ( i==MATROWS(omat) && j==MATCOLS(omat) )
+        {
+          Print("  0");
+        }
+        else
+        {
+          Print("  0, ");
+        }
+      }
+    }
 …
+  }
   Print(");\n");
+}
+}
 #endif
 //<-
 …
   int i;
   points = (onePointP *)Alloc( (count+1) * sizeof(onePointP) );
+  for ( i= 0; i <= max; i++ ) {
+  for ( i= 0; i <= max; i++ )
+  {
     points[i]= (onePointP)Alloc( sizeof(onePoint) );
     points[i]->point= (Coord_t *)Alloc0( (dim+2) * sizeof(Coord_t) );
 …
   int i;
   int fdim= lifted ? dim+1 : dim+2;
+  for ( i= 0; i <= max; i++ ) {
+  for ( i= 0; i <= max; i++ )
+  {
     Free( (ADDRESS) points[i]->point, fdim * sizeof(Coord_t) );
     Free( (ADDRESS) points[i], sizeof(onePoint) );
 …
 inline bool pointSet::checkMem()
+{
+  if ( num >= max ) {
+  if ( num >= max )
+  {
     int i;
     int fdim= lifted ? dim+1 : dim+2;
+    points= (onePointP*)ReAlloc( points,
+                                 (max+1) * sizeof(onePointP),
+                                 (2*max + 1) * sizeof(onePointP) );
+    for ( i= max+1; i <= max*2; i++ ) {
+    points= (onePointP*)ReAlloc( points,
+                                 (max+1) * sizeof(onePointP),
+                                 (2*max + 1) * sizeof(onePointP) );
+    for ( i= max+1; i <= max*2; i++ )
+    {
       points[i]= (onePointP)Alloc( sizeof(struct onePoint) );
       points[i]->point= (Coord_t *)Alloc0( fdim * sizeof(Coord_t) );
 …
+{
   assume( indx > 0 && indx <= num );
+  if ( indx != num ) {
+  if ( indx != num )
+  {
     onePointP tmp;
     tmp= points[indx];
 …
   int i,j;
+  for ( i= 1; i <= num; i++ ) {
+  for ( i= 1; i <= num; i++ )
+  {
     for ( j= 1; j <= dim; j++ )
       if ( points[i]->point[j] != vert->point[j] ) break;
     if ( j > dim ) break;
+  }
+  if ( i > num ) {
+  if ( i > num )
+  {
     addPoint( vert );
     return true;
 …
   int i,j;
+  for ( i= 1; i <= num; i++ ) {
+  for ( i= 1; i <= num; i++ )
+  {
     for ( j= 1; j <= dim; j++ )
       if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
     if ( j > dim ) break;
+  }
+  if ( i > num ) {
+  if ( i > num )
+  {
     addPoint( vert );
     return true;
 …
   vert= (Exponent_t *)Alloc( (dim+1) * sizeof(Exponent_t) );
+  while ( piter ) {
+  while ( piter )
+  {
     pGetExpV( piter, vert );
+    for ( i= 1; i <= num; i++ ) {
+    for ( i= 1; i <= num; i++ )
+    {
       for ( j= 1; j <= dim; j++ )
         if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
+        if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
       if ( j > dim ) break;
+    }
+    if ( i > num ) {
+    if ( i > num )
+    {
       addPoint( vert );
+    }
 …
   pGetExpV( p, vert );
+  for ( i= 1; i <= num; i++ ) {
+  for ( i= 1; i <= num; i++ )
+  {
     for ( j= 1; j <= dim; j++ )
       if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
 …
+{
   int i;
+  for ( i= 1; i <= dim; i++ ) {
+    if ( points[a]->point[i] > points[b]->point[i] ) {
+  for ( i= 1; i <= dim; i++ )
+  {
+    if ( points[a]->point[i] > points[b]->point[i] )
+    {
       return false;
+    }
+    if ( points[a]->point[i] < points[b]->point[i] ) {
+    if ( points[a]->point[i] < points[b]->point[i] )
+    {
       return true;
+    }
 …
+{
   int i;
+  for ( i= 1; i <= dim; i++ ) {
+    if ( points[a]->point[i] < points[b]->point[i] ) {
+  for ( i= 1; i <= dim; i++ )
+  {
+    if ( points[a]->point[i] < points[b]->point[i] )
+    {
       return false;
+    }
+    if ( points[a]->point[i] > points[b]->point[i] ) {
+    if ( points[a]->point[i] > points[b]->point[i] )
+    {
       return true;
+    }
 …
   onePointP tmp;
+  while ( found ) {
+  while ( found )
+  {
     found= false;
+    for ( i= 1; i < num; i++ ) {
+      if ( larger( i, i+1 ) ) {
+        tmp= points[i];
+        points[i]= points[i+1];
+        points[i+1]= tmp;
+        found= true;
+    for ( i= 1; i < num; i++ )
+    {
+      if ( larger( i, i+1 ) )
+      {
+        tmp= points[i];
+        points[i]= points[i+1];
+        points[i+1]= tmp;
+        found= true;
+      }
+    }
 …
   time_t *tp = NULL;
   seed = ((int)time(tp) % MAXSEED) + index; // secs since 1/1/70 ->[1,MAXSEED]
+  seed = ((int)time(tp) % MAXSEED) + index; // secs since 1/1/70 ->[1,MAXSEED]
   //seed= 100+index;
   dim++;
+  if ( !l ) {
+  if ( l==NULL )
+  {
     outerL= false;
     l= (int *)Alloc( (dim+1) * sizeof(int) ); // [1..dim-1]
     srand48( (long)(seed*(index+MAXSEED*18)) );  // index+MAXSEED*MAXVARS
+    for(i = 1; i < dim; i++) {
+    for(i = 1; i < dim; i++)
+    {
       l[i]= 1 + ((unsigned int) mrand48()) % LIFT_COOR;
+    }
+  }
+  for ( j=1; j <= num; j++ ) {
+  for ( j=1; j <= num; j++ )
+  {
     sum= 0;
+    for ( i=1; i < dim; i++ ) {
+    for ( i=1; i < dim; i++ )
+    {
       sum += (int)points[j]->point[i] * l[i];
+    }
     points[j]->point[dim]= sum;
+  }
 #ifdef mprDEBUG_ALL
   Print(" lift vector: ");
 …
 #ifdef mprDEBUG_ALL
   Print(" lifted points: \n");
+  for ( j=1; j <= num; j++ ) {
+  for ( j=1; j <= num; j++ )
+  {
     Print("%d: <",j);print_exp(points[j],dim);Print(">\n");
+  }
 …
+{
   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++) {
+  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;
+  }
   for( i= 1; i <= n; i++) {     // each row constraints one coor
+  for( i= 1; i <= n; i++) {        // each row constraints one coor
     pLP->LiPM[i+2][1] = (mprfloat)pGetExp(pointPoly,i);
     col = 2;
+    for( j= 1; j <= m; j++ ) {
+      if( j != site ) {
+    for( j= 1; j <= m; j++ )
+    {
+      if( j != site )
+      {
         pLP->LiPM[i+2][col] = -(mprfloat)pGetExp( monomAt(p,j), i );
         col++;
 …
 #endif
 #if 1
+  if ( numcons + 1 > pLP->LiPM_rows ) Werror("convexHull::inHull: #rows > #pLP->LiPM_rows!");
+  if ( numcols + 1 > pLP->LiPM_cols ) Werror("convexHull::inHull: #cols > #pLP->LiPM_cols!");
+  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);
+  return (icase == 0);
+}
 …
   Exponent_t * vert;
   n= pVariables;
+  n= pVariables;
   vert= (Exponent_t *)Alloc( (idelem+1) * sizeof(Exponent_t) );
   Q = (pointSet **)Alloc( idelem * sizeof(pointSet*) ); // support hulls
+  Q = (pointSet **)Alloc( idelem * sizeof(pointSet*) );        // support hulls
   for ( i= 0; i < idelem; i++ )
     Q[i] = new pointSet( pVariables, i+1, pLength((gls->m)[i])+1 );
+  for( i= 0; i < idelem; i++ ) {
+  for( i= 0; i < idelem; i++ )
+  {
     k=1;
     m = pLength( (gls->m)[i] );
     poly p= (gls->m)[i];
+    for( j= 1; j <= m; j++) {  // für jeden Exponentvektor
+      if( !inHull( (gls->m)[i], p, m, j ) ) {
+        pGetExpV( p, vert );
+    for( j= 1; j <= m; j++) {  // für jeden Exponentvektor
+      if( !inHull( (gls->m)[i], p, m, j ) )
+      {
+        pGetExpV( p, vert );
         Q[i]->addPoint( vert );
+        k++;
+        mprSTICKYPROT(ST_SPARSE_VADD);
+      } else {
+        mprSTICKYPROT(ST_SPARSE_VREJ);
+        k++;
+        mprSTICKYPROT(ST_SPARSE_VADD);
+      }
+      else
+      {
+        mprSTICKYPROT(ST_SPARSE_VREJ);
+      }
       pIter( p );
 …
 #ifdef mprDEBUG_PROT
   PrintLn();
+  for( i= 0; i < idelem; i++ ) {
+  for( i= 0; i < idelem; i++ )
+  {
     Print(" \\Conv(Qi[%d]): #%d\n", i,Q[i]->num );
+    for ( j=1; j <= Q[i]->num; j++ ) {
+    for ( j=1; j <= Q[i]->num; j++ )
+    {
       Print("%d: <",j);print_exp( (*Q[i])[j] , pVariables );Print(">\n");
+    }
 …
   numverts = 0;
+  for( i=0; i<=n; i++) {
+  for( i=0; i<=n; i++)
+  {
     numverts += Qi[i]->num;
+  }
+  cols = numverts + 2;
+  //if( dim < 1 || dim > n ) Werror("mayanPyramidAlg::vDistance: Known coords dim off range");
+  pLP->LiPM[1][1] = 0.0;
+  pLP->LiPM[1][2] = 1.0;        // maximize
+  cols = numverts + 2;
+  //if( dim < 1 || dim > n )
+  //  WerrorS("mayanPyramidAlg::vDistance: Known coords dim off range");
+  pLP->LiPM[1][1] = 0.0;
+  pLP->LiPM[1][2] = 1.0;        // maximize
   for( j=3; j<=cols; j++) pLP->LiPM[1][j] = 0.0;
   for( i=0; i <= n; i++ ) {
+  for( i=0; i <= n; i++ ) {
     pLP->LiPM[i+2][1] = 1.0;
     pLP->LiPM[i+2][2] = 0.0;
+  }
   for( i=1; i<=dim; i++) {
+  for( i=1; i<=dim; i++) {
     pLP->LiPM[n+2+i][1] = (mprfloat)(acoords[i-1]);
     pLP->LiPM[n+2+i][2] = -shift[i];
 …
   ii = -1;
   col = 2;
   for ( i= 0; i <= n; i++ ) {
+  for ( i= 0; i <= n; i++ ) {
     ii++;
+    for( k= 1; k <= Qi[ii]->num; k++ ) {
+    for( k= 1; k <= Qi[ii]->num; k++ )
+    {
       col++;
+      for ( r= 0; r <= n; r++ ) {
+        if ( r == i ) pLP->LiPM[r+2][col] = -1.0;
+        else pLP->LiPM[r+2][col] = 0.0;
+      for ( r= 0; r <= n; r++ )
+      {
+        if ( r == i ) pLP->LiPM[r+2][col] = -1.0;
+        else pLP->LiPM[r+2][col] = 0.0;
+      }
       for( r= 1; r <= dim; r++ ) pLP->LiPM[r+n+2][col] = -(mprfloat)((*Qi[ii])[k]->point[r]);
 …
+  }
+  if( col != cols) Werror("mayanPyramidAlg::vDistance:"
+                          "setting up matrix for udist: col %d != cols %d",col,cols);
+  constr = n+dim+1;
+  if( col != cols)
+    Werror("mayanPyramidAlg::vDistance:"
+           "setting up matrix for udist: col %d != cols %d",col,cols);
+  constr = n+dim+1;
 #ifdef mprDEBUG_ALL
 …
 #if 1
+  if ( constr + 1 > pLP->LiPM_rows ) Werror("mayanPyramidAlg::vDistance: #rows > #pLP->LiPM_rows!");
+  if ( cols + 1 > pLP->LiPM_cols ) Werror("mayanPyramidAlg::vDistance: #cols > #pLP->LiPM_cols!");
+  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
 …
   print_bmat( pLP->LiPM, constr+1, cols+1-constr, cols, pLP->iposv);
 #endif
   if( icase != 0 ) {  // check for errors
+    Werror("mayanPyramidAlg::vDistance: \n");
+    if( icase == 1 ) Werror(" vDistance: Unbounded v-distance: probably 1st v-coor=0");
+    if( icase == -1 ) Werror(" vDistance: Infeasible v-distance");
+    WerrorS("mayanPyramidAlg::vDistance:");
+    if( icase == 1 )
+      WerrorS(" vDistance: Unbounded v-distance: probably 1st v-coor=0");
+    if( icase == -1 )
+      WerrorS(" vDistance: Infeasible v-distance");
     return -1.0;
+  }
 …
   // common part of the matrix
   pLP->LiPM[1][1] = 0.0;
   for( i=2; i<=n+2; i++) {
     pLP->LiPM[i][1] = 1.0;      // 1st col
     pLP->LiPM[i][2] = 0.0;      // 2nd col
+  pLP->LiPM[1][1] = 0.0;
+  for( i=2; i<=n+2; i++) {
+    pLP->LiPM[i][1] = 1.0;        // 1st col
+    pLP->LiPM[i][2] = 0.0;        // 2nd col
+  }
   la_cons_row = 1;
   cols = 2;
+  for( i=0; i<=n; i++) {
+  for( i=0; i<=n; i++)
+  {
     la_cons_row++;
+    for( j=1; j<= Qi[i]->num; j++) {
+    for( j=1; j<= Qi[i]->num; j++)
+    {
       cols++;
       pLP->LiPM[1][cols] = 0.0; // set 1st row 0
+      pLP->LiPM[1][cols] = 0.0;        // set 1st row 0
       for( k=2; k<=n+2; k++) {  // lambdas sum up to 1
         if( k != la_cons_row) pLP->LiPM[k][cols] = 0.0;
         else pLP->LiPM[k][cols] = -1.0;
+        if( k != la_cons_row) pLP->LiPM[k][cols] = 0.0;
+        else pLP->LiPM[k][cols] = -1.0;
+      }
       for( k=1; k<=n; k++) pLP->LiPM[k+n+2][cols] = -(mprfloat)((*Qi[i])[j]->point[k]);
 …
   } // i
   for( i= 0; i < dim; i++ ) {           // fixed coords
+  for( i= 0; i < dim; i++ ) {                // fixed coords
     pLP->LiPM[i+n+3][1] = acoords[i];
     pLP->LiPM[i+n+3][2] = 0.0;
 …
   pLP->LiPM[dim+n+3][1] = 0.0;
   pLP->LiPM[1][2] = -1.0;                       // minimize
+  pLP->LiPM[1][2] = -1.0;                        // minimize
   pLP->LiPM[dim+n+3][2] = 1.0;
 …
 #if 1
+  if ( cons + 1 > pLP->LiPM_rows ) Werror(" mn_mx_MinkowskiSum: #rows > #pLP->LiPM_rows!");
+  if ( cols + 1 > pLP->LiPM_cols ) Werror(" mn_mx_MinkowskiSum: #cols > #pLP->LiPM_cols!");
+  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) Werror(" mn_mx_MinkowskiSum: LinearProgram: minR: infeasible");
+    if( icase > 0) Werror(" mn_mx_MinkowskiSum: LinearProgram: minR: unbounded");
+    if( icase < 0)
+      WerrorS(" mn_mx_MinkowskiSum: LinearProgram: minR: infeasible");
+    else if( icase > 0)
+      WerrorS(" mn_mx_MinkowskiSum: LinearProgram: minR: unbounded");
+  }
 …
   // common part of the matrix again
   pLP->LiPM[1][1] = 0.0;
   for( i=2; i<=n+2; i++) {
     pLP->LiPM[i][1] = 1.0;
     pLP->LiPM[i][2] = 0.0;
+  pLP->LiPM[1][1] = 0.0;
+  for( i=2; i<=n+2; i++) {
+    pLP->LiPM[i][1] = 1.0;
+    pLP->LiPM[i][2] = 0.0;
+  }
   la_cons_row = 1;
   cols = 2;
+  for( i=0; i<=n; i++) {
+  for( i=0; i<=n; i++)
+  {
     la_cons_row++;
+    for( j=1; j<=Qi[i]->num; j++) {
+    for( j=1; j<=Qi[i]->num; j++)
+    {
       cols++;
       pLP->LiPM[1][cols] = 0.0;
       for( k=2; k<=n+2; k++) {
         if( k != la_cons_row) pLP->LiPM[k][cols] = 0.0;
         else pLP->LiPM[k][cols] = -1.0;
+      pLP->LiPM[1][cols] = 0.0;
+      for( k=2; k<=n+2; k++) {
+        if( k != la_cons_row) pLP->LiPM[k][cols] = 0.0;
+        else pLP->LiPM[k][cols] = -1.0;
+      }
       for( k=1; k<=n; k++) pLP->LiPM[k+n+2][cols] = -(mprfloat)((*Qi[i])[j]->point[k]);
 …
   }  // i
   for( i= 0; i < dim; i++ ) {           // fixed coords
+  for( i= 0; i < dim; i++ ) {                // fixed coords
     pLP->LiPM[i+n+3][1] = acoords[i];
     pLP->LiPM[i+n+3][2] = 0.0;
 …
   pLP->LiPM[dim+n+3][1] = 0.0;
   pLP->LiPM[1][2] = 1.0;                        // maximize
   pLP->LiPM[dim+n+3][2] = 1.0;          // var = sum of pnt coords
+  pLP->LiPM[1][2] = 1.0;                        // maximize
+  pLP->LiPM[dim+n+3][2] = 1.0;                // var = sum of pnt coords
 #ifdef mprDEBUG_ALL
 …
 #if 1
+  if ( cons + 1 > pLP->LiPM_rows ) Werror(" mn_mx_MinkowskiSum: #rows > #pLP->LiPM_rows!");
+  if ( cols + 1 > pLP->LiPM_cols ) Werror(" mn_mx_MinkowskiSum: #cols > #pLP->LiPM_cols!");
+  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(pLP->LiPM, cons, cols-1, 0, 0, cons, &icase, pLP->izrov, pLP->iposv);
+  if ( icase != 0 ) {
+    if( icase < 0) Werror(" mn_mx_MinkowskiSum: LinearProgram: maxR: infeasible");
+    if( icase > 0) Werror(" mn_mx_MinkowskiSum: LinearProgram: maxR: unbounded");
+  }
+  *maxR = (Coord_t)( pLP->LiPM[1][1] + SIMPLEX_EPS );
+  if ( icase != 0 )
+  {
+    if( icase < 0)
+      WerrorS(" mn_mx_MinkowskiSum: LinearProgram: maxR: infeasible");
+    else if( icase > 0)
+      WerrorS(" mn_mx_MinkowskiSum: LinearProgram: maxR: unbounded");
+  }
+  *maxR = (Coord_t)( pLP->LiPM[1][1] + SIMPLEX_EPS );
 #ifdef mprDEBUG_ALL
 …
 // mprSTICKYPROT:
 // ST_SPARSE_VREJ: rejected point
 // ST_SPARSE_VADD: added point to set
+// ST_SPARSE_VREJ: rejected point
+// ST_SPARSE_VADD: added point to set
 bool mayanPyramidAlg::storeMinkowskiSumPoint()
+{
 …
   dist= vDistance( &(acoords[0]), n );
+  // store only points with v-distance > minVdist
+  if( dist <= MINVDIST + SIMPLEX_EPS ) {
+  // store only points with v-distance > minVdist
+  if( dist <= MINVDIST + SIMPLEX_EPS )
+  {
     mprSTICKYPROT(ST_SPARSE_VREJ);
     return false;
 …
   E->addPoint( &(acoords[0]) );
   mprSTICKYPROT(ST_SPARSE_VADD);
   return true;
+  return true;
+}
 …
   // step 5 -> terminate
   if( dim == n-1 ) {
     int lastKilled = 0;
+  if( dim == n-1 ) {
+    int lastKilled = 0;
     // insert points
     acoords[dim] = minR;
+    while( acoords[dim] <= maxR ) {
+    while( acoords[dim] <= maxR )
+    {
       if( !storeMinkowskiSumPoint() ) lastKilled++;
       acoords[dim]++;
 …
   // step 4 -> recurse at step 3
+  acoords[dim] = minR;
+  while ( acoords[dim] <= maxR ) {
+  acoords[dim] = minR;
+  while ( acoords[dim] <= maxR )
+  {
     if ( (acoords[dim] > minR) && (acoords[dim] <= maxR) ) {     // acoords[dim] >= minR  ??
       mprSTICKYPROT(ST_SPARSE_MREC1);
       runMayanPyramid( dim + 1 );         // recurse with higer dimension
     } else {
+    } else {
       // get v-distance of pt
       dist= vDistance( &(acoords[0]), dim + 1 );// dim+1 == known coordinates
+      if( dist >= SIMPLEX_EPS ) {
+        mprSTICKYPROT(ST_SPARSE_MREC2);
+        runMayanPyramid( dim + 1 );       // recurse with higer dimension
+      if( dist >= SIMPLEX_EPS )
+      {
+        mprSTICKYPROT(ST_SPARSE_MREC2);
+        runMayanPyramid( dim + 1 );       // recurse with higer dimension
+      }
+    }
 …
   int i,n= pVariables;
   int loffset= 0;
+  for ( i= 0; i <= n; i++ ) {
+    if ( (loffset < indx) && (indx <= pQ[i]->num + loffset) ) {
+  for ( i= 0; i <= n; i++ )
+  {
+    if ( (loffset < indx) && (indx <= pQ[i]->num + loffset) )
+    {
       *set= i;
       *pnt= indx-loffset;
 …
   // fill in LP matrix
+  for ( i= 0; i <= n; i++ ) {
+  for ( i= 0; i <= n; i++ )
+  {
     size= pQ[i]->num;
+    for ( k= 1; k <= size; k++ ) {
+    for ( k= 1; k <= size; k++ )
+    {
       ncols++;
       // objective funtion, minimize
       LP.LiPM[1][ncols] = - ( (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;
         else
           LP.LiPM[j+2][ncols] = 0.0;
+      for ( j = 0; j <= n; j++ )
+        if ( i==j )
+          LP.LiPM[j+2][ncols] = -1.0;
+        else
+          LP.LiPM[j+2][ncols] = 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;
+  for ( j= 1; j <= n; j++ ) {
+      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;
+  for ( j= 1; j <= n; j++ )
+  {
     LP.LiPM[j+n+2][1]= (mprfloat)(*E)[vert]->point[j] - shift[j];
+  }
 …
   PrintLn();
   Print(" n= %d, NumCons=M= %d, ncols=N= %d\n",n,NumCons,ncols);
   print_mat(LP.LiPM, NumCons+1, ncols+1);
 #endif
+  print_mat(LP.LiPM, NumCons+1, ncols+1);
+#endif
 #if 1
+  if ( NumCons + 1 > LP.LiPM_rows ) Werror("resMatrixSparse::RC: #rows > #LP.LiPM_rows!");
+  if ( ncols + 1 > LP.LiPM_cols ) Werror("resMatrixSparse::RC: #cols > #LP.LiPM_cols!");
+  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
 …
   simplx(LP.LiPM, NumCons, ncols,  0,  0, NumCons, &icase, LP.izrov, LP.iposv);
+  if (icase < 0) {
+  if (icase < 0)
+  {
     // infeasibility: the point does not lie in a cell -> remove it
     return -1;
 …
   //print_mat(LP.LiPM, NumCons+1, ncols);
 #endif
 #if 1
   // sort LP results
+  while (found) {
+  while (found)
+  {
     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;
+        found= true;
+    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;
+        found= true;
+      }
+    }
 …
   c=0;
   optSum= (setID*)Alloc( (NumCons) * sizeof(struct setID) );
+  for ( i= 0; i < NumCons; i++ ) {
+  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]);
+        Werror(" resMatrixSparse::RC: remapXiToPoint faild!");
+        return -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;
+      }
       bucket[optSum[c].set]++;
 …
   // find last min in bucket[]: maximum i such that Fi is a point
   c= 0;
+  for ( i= 1; i < E->dim; i++ ) {
+    if ( bucket[c] >= bucket[i] ) {
+  for ( i= 1; i < E->dim; i++ )
+  {
+    if ( bucket[c] >= bucket[i] )
+    {
       c= i;
+    }
+  }
   // find matching point set
+  for ( i= onum - 1; i >= 0; i-- ) {
+    if ( optSum[i].set == c )
+  for ( i= onum - 1; i >= 0; i-- )
+  {
+    if ( optSum[i].set == c )
       break;
+  }
 …
   (*E)[vert]->rc.pnt= optSum[i].pnt;
   (*E)[vert]->rcPnt= (*pQ[c])[optSum[i].pnt];
   // count
+  // count
   if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
 #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++ ) {
+  for ( j= 0; j < E->dim; j++ )
+  {
     Print(" %d",bucket[j]);
+  }
   Print(" }\n optimal Sum: Qi ");
+  for ( j= 0; j < NumCons; j++ ) {
+  for ( j= 0; j < NumCons; j++ )
+  {
     Print(" [ %d, %d ]",optSum[j].set,optSum[j].pnt);
+  }
 …
   mprSTICKYPROT(ST_SPARSE_RC);
   return (int)(-LP.LiPM[1][1] * SCALEDOWN);
+  return (int)(-LP.LiPM[1][1] * SCALEDOWN);
+}
 …
   for ( i= 1; i <= E->num; i++ ) {       // for every row
     E->getRowMP( i, epp_mon );           // compute (p-a[ij]), (i,j) = RC(p)
     pSetExpV( epp, epp_mon );
     //
+    pSetExpV( epp, epp_mon );
+    //
     rowp= pMult( pCopy(epp), pCopy((gls->m)[(*E)[i]->rc.set]) );  // x^(p-a[ij]) * f(i)
     pSetm( rowp );
 …
     // get column for every monomial in rowp and store it
     iterp= rowp;
+    while ( iterp ) {
+    while ( iterp )
+    {
       epos= E->getExpPos( iterp );
+      if ( epos == 0 ) {
+        // this can happen, if the shift vektor or the lift funktions
+        // are not generically choosen.
+        Werror("resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
+               i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
+        return i;
+      if ( epos == 0 )
+      {
+        // this can happen, if the shift vektor or the lift funktions
+        // are not generically choosen.
+        Werror("resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
+               i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
+        return i;
+      }
       pSetExpV(iterp,eexp);
       pSetComp(iterp, epos );
       if ( (*E)[i]->rc.set == linPolyS ) { // store coeff positions
         IMATELEM(*uRPos,rp,cp)= epos;
         cp++;
+        IMATELEM(*uRPos,rp,cp)= epos;
+        cp++;
+      }
       pIter( iterp );
 …
 #ifdef mprDEBUG_ALL
+  if ( E->num <= 40 ) {
+  if ( E->num <= 40 )
+  {
     matrix mout= idModule2Matrix( idCopy(rmat) );
     print_matrix(mout);
+  }
+  for ( i= 1; i <= numSet0; i++ ) {
+  for ( i= 1; i <= numSet0; i++ )
+  {
     Print(" row  %d contains coeffs of f_%d\n",IMATELEM(*uRPos,i,1),linPolyS);
+  }
 …
   srand((long)(((int)time(tp) % MAXSEED)+MAXSEED*18));
+  while ( i <= dim ) {
+  while ( i <= dim )
+  {
     shift[i]= (mprfloat) (RVMULT*rand()/(RAND_MAX+1.0));
     i++;
+    for ( j= 1; j < i-1; j++ ) {
+      if ( (shift[j] < shift[i-1] + SIMPLEX_EPS) && (shift[j] > shift[i-1] - SIMPLEX_EPS) ) {
+        i--;
+        break;
+    for ( j= 1; j < i-1; j++ )
+    {
+      if ( (shift[j] < shift[i-1] + SIMPLEX_EPS) && (shift[j] > shift[i-1] - SIMPLEX_EPS) )
+      {
+        i--;
+        break;
+      }
+    }
 …
   vs= new pointSet( dim );
+  for ( j= 1; j <= Q1->num; j++ ) {
+    for ( k= 1; k <= Q2->num; k++ ) {
+      for ( l= 1; l <= dim; l++ ) {
+        vert.point[l]= (*Q1)[j]->point[l] + (*Q2)[k]->point[l];
+  for ( j= 1; j <= Q1->num; j++ )
+  {
+    for ( k= 1; k <= Q2->num; k++ )
+    {
+      for ( l= 1; l <= dim; l++ )
+      {
+        vert.point[l]= (*Q1)[j]->point[l] + (*Q2)[k]->point[l];
+      }
       vs->mergeWithExp( &vert );
 …
   for ( j= 1; j <= pQ[0]->num; j++ ) vs->addPoint( (*pQ[0])[j] );
+  for ( j= 1; j < numq; j++ ) {
+  for ( j= 1; j < numq; j++ )
+  {
     vs_old= vs;
     vs= minkSumTwo( vs_old, pQ[j], dim );
 …
   mprfloat shift[MAXVARS+2];   // shiftvector delta, index [1..dim]
+  if ( pVariables > MAXVARS ) {
+    Werror("resMatrixSparse::resMatrixSparse: To many variables!");
+  if ( pVariables > MAXVARS )
+  {
+    WerrorS("resMatrixSparse::resMatrixSparse: To many variables!");
     return;
+  }
 …
   // 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++ ) {
+  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) );
 …
 #ifdef mprMINKSUM
   E= minkSumAll( Qi, n+1, n);
+  E= minkSumAll( Qi, n+1, n);
 #else
   // get inner points
 …
 #endif
   Print("\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
+  for ( pnt= 1; pnt <= E->num; pnt++ ) {
+  for ( pnt= 1; pnt <= E->num; pnt++ )
+  {
     Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );Print(">\n");
+  }
 …
 #ifdef mprTEST
   int lift[5][5];
   lift[0][1]=3; lift[0][2]=4; lift[0][3]=8;  lift[0][4]=2;
   lift[1][1]=6; lift[1][2]=1; lift[1][3]=7;  lift[1][4]=4;
+  lift[0][1]=3; lift[0][2]=4; lift[0][3]=8;  lift[0][4]=2;
+  lift[1][1]=6; lift[1][2]=1; lift[1][3]=7;  lift[1][4]=4;
   lift[2][1]=2; lift[2][2]=5; lift[2][3]=9;  lift[2][4]=6;
   lift[3][1]=2; lift[3][2]=1; lift[3][3]=9;  lift[3][4]=5;
 …
   // run Row Content Function for every point in E
+  for ( pnt= 1; pnt <= E->num; pnt++ ) {
+  for ( pnt= 1; pnt <= E->num; pnt++ )
+  {
     RC( Qi, E, pnt, shift );
+  }
 …
   // remove points not in cells
   k= E->num;
+  for ( pnt= k; pnt > 0; pnt-- ) {
+    if ( (*E)[pnt]->rcPnt == NULL ) {
+  for ( pnt= k; pnt > 0; pnt-- )
+  {
+    if ( (*E)[pnt]->rcPnt == NULL )
+    {
       E->removePoint(pnt);
       mprSTICKYPROT(ST_SPARSE_RCRJ);
 …
 #ifdef mprDEBUG_PROT
   Print(" points which lie in a cell:\n");
+  for ( pnt= 1; pnt <= E->num; pnt++ ) {
+  for ( pnt= 1; pnt <= E->num; pnt++ )
+  {
     Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );Print(">\n");
+  }
 …
 #ifdef mprDEBUG_PROT
   Print(" points with a[ij] (%d):\n",E->num);
+  for ( pnt= 1; pnt <= E->num; pnt++ ) {
+  for ( pnt= 1; pnt <= E->num; pnt++ )
+  {
     Print("Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->dim );
     Print(">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
 …
   // now create matrix
+  if ( createMatrix( E ) != E->num ) {
+  if ( createMatrix( E ) != E->num )
+  {
     // this can happen if the shiftvector shift is to large or not generic
     istate= resMatrixBase::fatalError;
     Werror("resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
+    WerrorS("resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
     goto theEnd;
     //return;
 …
  theEnd:
   // clean up
+  for ( i= 0; i < idelem; i++ ) {
+  for ( i= 0; i < idelem; i++ )
+  {
     delete Qi[i];
+  }
 …
   delete E;
+  for( i= 0; i < LP.LiPM_rows; i++ ) {
+  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.LiPM, LP.LiPM_rows * sizeof(mprfloat *) );
   Free( (ADDRESS) LP.iposv, (idelem * MAXPOINTS) * sizeof(int) );
 …
   // now fill in coeffs of f0
+  for ( i= 1; i <= numSet0; i++ ) {
+  for ( i= 1; i <= numSet0; i++ )
+  {
     pgls= (gls->m)[0]; // f0
 …
     // u_1,..,u_k
     cp=2;
+    while ( pNext(pgls) ) {
+    while ( pNext(pgls) )
+    {
       phelp= pOne();
       pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
       pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
+      if ( piter ) {
+        pNext(piter)= phelp;
+        piter= phelp;
+      } else {
+        pp= phelp;
+        piter= phelp;
+      if ( piter )
+      {
+        pNext(piter)= phelp;
+        piter= phelp;
+      }
+      else
+      {
+        pp= phelp;
+        piter= phelp;
+      }
       cp++;
 …
   mprPROTnl("smCallDet");
+  for ( i= 1; i <= numSet0; i++ ) {
+  for ( i= 1; i <= numSet0; i++ )
+  {
     pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
     pDelete( &pp );
 …
     piter= NULL;
     // u_1,..,u_n
+    for ( cp= 2; cp <= idelem; cp++ ) {
+      if ( !nIsZero(evpoint[cp-1]) ) {
+        phelp= pOne();
+        pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
+        pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
+        if ( piter ) {
+          pNext(piter)= phelp;
+          piter= phelp;
+        } else {
+          pp= phelp;
+          piter= phelp;
+        }
+    for ( cp= 2; cp <= idelem; cp++ )
+    {
+      if ( !nIsZero(evpoint[cp-1]) )
+      {
+        phelp= pOne();
+        pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
+        pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
+        if ( piter )
+        {
+          pNext(piter)= phelp;
+          piter= phelp;
+        }
+        else
+        {
+          pp= phelp;
+          piter= phelp;
+        }
+      }
+    }
 …
   mprPROTnl("smCallDet");
+  for ( i= 1; i <= numSet0; i++ ) {
+  for ( i= 1; i <= numSet0; i++ )
+  {
     pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
     pDelete( &pp );
 …
     piter= NULL;
     for ( cp= 2; cp <= idelem; cp++ ) { // u1 .. un
+      if ( !nIsZero(evpoint[cp-1]) ) {
+        phelp= pOne();
+        pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
+        pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
+        if ( piter ) {
+          pNext(piter)= phelp;
+          piter= phelp;
+        } else {
+          pp= phelp;
+          piter= phelp;
+        }
+      if ( !nIsZero(evpoint[cp-1]) )
+      {
+        phelp= pOne();
+        pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
+        pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
+        if ( piter )
+        {
+          pNext(piter)= phelp;
+          piter= phelp;
+        }
+        else
+        {
+          pp= phelp;
+          piter= phelp;
+        }
+      }
+    }
 …
    * _gls: system of multivariate polynoms
    * special: -1 -> resMatrixDense is a symbolic matrix
    *    0,1, ... -> resMatrixDense ist eine u-Resultante, wobei special das
+   *    0,1, ... -> resMatrixDense ist eine u-Resultante, wobei special das
    *                        lineare u-Polynom angibt
    */
 …
    */
   const number getSubDet();
 private:
   /** deactivated copy constructor */
 …
   void generateMonomData( int deg, intvec* polyDegs , intvec* iVO );
   /** Recursively generate all homogeneous monoms of
+  /** Recursively generate all homogeneous monoms of
    * pVariables variables of degree deg.
    */
 …
 //<-
 //-> struct resVector
+//-> struct resVector
 /* Holds a row vector of the dense resultant matrix */
 struct resVector
+{
 public:
+  void init() {
+  void init()
+  {
     isReduced = FALSE;
     elementOfS = SFREE;
+    elementOfS = SFREE;
     mon = NULL;
+  }
+  void init( const poly m ) {
+  void init( const poly m )
+  {
     isReduced = FALSE;
     elementOfS = SFREE;
     mon = m;
+    elementOfS = SFREE;
+    mon = m;
+  }
   /** index von 0 ... numVectors-1 */
   poly getElem( const int i );
+  poly getElem( const int i );
   /** index von 0 ... numVectors-1 */
 …
   int elementOfS;
   /** holds the index of u0, u1, ..., un, if (elementOfS == linPolyS)
    *  the size is given by pVariables
+  /** holds the index of u0, u1, ..., un, if (elementOfS == linPolyS)
+   *  the size is given by pVariables
    */
   int *numColParNr;
+  int *numColParNr;
   /** holds the column vector if (elementOfS != linPolyS) */
   number *numColVector;
+  number *numColVector;
   /** size of numColVector */
   int numColVectorSize;
+  int numColVectorSize;
   number *numColVecCopy;
 …
 //-> resVector::*
 poly resVector::getElem( const int i ) // inline ???
+{
+{
   assume( 0 < i || i > numColVectorSize );
   poly out= pOne();
   pSetCoeff( out, numColVector[i] );
   pTest( out );
   return( out );
+  return( out );
+}
 …
+{
   assume( i >= 0 && i < numColVectorSize );
   return( numColVector[i] );
+  return( numColVector[i] );
+}
 //<-
 …
 resMatrixDense::resMatrixDense( const ideal _gls, const int special )
   : resMatrixBase()
+{
+{
   int i;
   sourceRing=currRing;
   gls= idCopy( _gls );
 …
   // init all
   generateBaseData();
+  generateBaseData();
   totDeg= 1;
+  for ( i= 0; i < IDELEMS(gls); i++ ) {
+  for ( i= 0; i < IDELEMS(gls); i++ )
+  {
     totDeg*=pTotaldegree( (gls->m)[i] );
+  }
 …
   istate= resMatrixBase::ready;
+}
 resMatrixDense::~resMatrixDense()
+{
   int i,j;
+  for (i=0; i < numVectors; i++)
+    {
+      pDelete( &resVectorList[i].mon );
+      pDelete( &resVectorList[i].dividedBy );
+      for ( j=0; j < resVectorList[i].numColVectorSize; j++ ) {
+        nDelete( resVectorList[i].numColVector+j );
+      }
+      Free( (ADDRESS)resVectorList[i].numColVector, numVectors*sizeof( number ) );
+      Free( (ADDRESS)resVectorList[i].numColParNr, (pVariables+1) * sizeof(int) );
+    }
+  for (i=0; i < numVectors; i++)
+  {
+    pDelete( &resVectorList[i].mon );
+    pDelete( &resVectorList[i].dividedBy );
+    for ( j=0; j < resVectorList[i].numColVectorSize; j++ )
+    {
+        nDelete( resVectorList[i].numColVector+j );
+    }
+    Free( (ADDRESS)resVectorList[i].numColVector, numVectors*sizeof( number ) );
+    Free( (ADDRESS)resVectorList[i].numColParNr, (pVariables+1) * sizeof(int) );
+  }
   Free( (ADDRESS)resVectorList, veclistmax*sizeof( resVector ) );
+  // free matrix m
+  if ( m != NULL ) {
+  // free matrix m
+  if ( m != NULL )
+  {
     for ( i= 1; i <= numVectors; i++ )
       for ( j= 1; j <= numVectors; j++ )
         pDelete( &MATELEM(m , i, j) );
+        pDelete( &MATELEM(m , i, j) );
     Free( (ADDRESS)m->m, numVectors * numVectors * sizeof(poly) );
     Free( (ADDRESS)m, sizeof(ip_smatrix) );
 …
   int k,i,j;
   resVector *vecp;
   m= mpNew( numVectors, numVectors );
+  for ( i= 1; i <= MATROWS( m ); i++ )
+    for ( j= 1; j <= MATCOLS( m ); j++ ) {
+  for ( i= 1; i <= MATROWS( m ); i++ )
+    for ( j= 1; j <= MATCOLS( m ); j++ )
+    {
       MATELEM(m,i,j)= pInit();
       pSetCoeff0( MATELEM(m,i,j), nInit(0) );
+    }
+  for ( k= 0; k <= numVectors - 1; k++ ) {
+    if ( linPolyS == getMVector(k)->elementOfS ) {
+  for ( k= 0; k <= numVectors - 1; k++ )
+  {
+    if ( linPolyS == getMVector(k)->elementOfS )
+    {
       mprSTICKYPROT(ST_DENSE_FR);
+      for ( i= 0; i < pVariables; i++ ) {
+        MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i])= pInit();
+      }
+    } else {
+      for ( i= 0; i < pVariables; i++ )
+      {
+        MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i])= pInit();
+      }
+    }
+    else
+    {
       mprSTICKYPROT(ST_DENSE_NR);
       vecp= getMVector(k);
+      for ( i= 0; i < numVectors; i++)
+        if ( !nIsZero( vecp->getElemNum(i) ) ) {
+          MATELEM(m,numVectors - k,i + 1)= pInit();
+          pSetCoeff( MATELEM(m,numVectors - k,i + 1), nCopy(vecp->getElemNum(i)) );
+        }
+      for ( i= 0; i < numVectors; i++)
+        if ( !nIsZero( vecp->getElemNum(i) ) )
+        {
+          MATELEM(m,numVectors - k,i + 1)= pInit();
+          pSetCoeff( MATELEM(m,numVectors - k,i + 1), nCopy(vecp->getElemNum(i)) );
+        }
+    }
   } // for
   mprSTICKYPROT("\n");
 #ifdef mprDEBUG_ALL
+  for ( k= numVectors - 1; k >= 0; k-- )
+    if ( linPolyS == getMVector(k)->elementOfS ) {
+      for ( i=0; i < pVariables; i++ ) {
+        Print(" %d ",(getMVector(k)->numColParNr)[i]);
+  for ( k= numVectors - 1; k >= 0; k-- )
+    if ( linPolyS == getMVector(k)->elementOfS )
+    {
+      for ( i=0; i < pVariables; i++ )
+      {
+        Print(" %d ",(getMVector(k)->numColParNr)[i]);
+      }
       PrintLn();
+    }
+  for (i=1; i <= numVectors; i++) {
+    for (j=1; j <= numVectors; j++ ) {
+  for (i=1; i <= numVectors; i++)
+  {
+    for (j=1; j <= numVectors; j++ )
+    {
       pWrite0(MATELEM(m,i,j));Print("  ");
+    }
 …
     poly mon = pCopy( m );
     pSetm( mon );
+    if ( numVectors == veclistmax ) {
+      resVectorList= (resVector * )ReAlloc( resVectorList,
+                                            (veclistmax) * sizeof( resVector ),
+                                            (veclistmax + veclistblock) * sizeof( resVector ) );
+    if ( numVectors == veclistmax )
+    {
+      resVectorList= (resVector * )ReAlloc( resVectorList,
+                                            (veclistmax) * sizeof( resVector ),
+                                            (veclistmax + veclistblock) * sizeof( resVector ) );
       int k;
       for ( k= veclistmax; k < (veclistmax + veclistblock); k++ )
         resVectorList[k].init();
+      for ( k= veclistmax; k < (veclistmax + veclistblock); k++ )
+        resVectorList[k].init();
       veclistmax+= veclistblock;
       mprSTICKYPROT(ST_DENSE_MEM);
 …
     resVectorList[numVectors].init( mon );
     numVectors++;
     mprSTICKYPROT(ST_DENSE_NMON);
+    mprSTICKYPROT(ST_DENSE_NMON);
     return;
+  } else {
+  }
+  else
+  {
     if ( var == pVariables+1 ) return;
     poly newm = pCopy( m );
+    while ( deg >= 0 ) {
+    while ( deg >= 0 )
+    {
       generateMonoms( newm, var+1, deg );
       pIncrExp( newm, var );
 …
     pDelete( & newm );
+  }
   return;
+}
 …
   veclistmax= veclistblock;
   resVectorList= (resVector *)Alloc( veclistmax*sizeof( resVector ) );
   // Init resVector()s
   for ( j= veclistmax - 1; j >= 0; j-- ) resVectorList[j].init();
   numVectors= 0;
   // Generate all monoms of degree deg
+  // Generate all monoms of degree deg
   poly start= pOne();
   generateMonoms( start, 1, deg );
 …
   mprSTICKYPROT("\n");
   // Check for reduced monoms
   // First generate polyDegs.rows() monoms
   //  x(k)^(polyDegs[k]),  0 <= k < polyDegs.rows()
   ideal pDegDiv= idInit( polyDegs->rows(), 1 );
+  for ( k= 0; k < polyDegs->rows(); k++ ) {
+  for ( k= 0; k < polyDegs->rows(); k++ )
+  {
     poly p= pOne();
     pSetExp( p, k + 1, (*polyDegs)[k] );
 …
   // exactly one x(k)^(polyDegs[k]) that divides the monom.
   int divCount;
+  for ( j= numVectors - 1; j >= 0; j-- ) {
+  for ( j= numVectors - 1; j >= 0; j-- )
+  {
     divCount= 0;
     for ( k= 0; k < IDELEMS(pDegDiv); k++ )
       if ( pDivisibleBy2( (pDegDiv->m)[k], resVectorList[j].mon ) )
         divCount++;
+    for ( k= 0; k < IDELEMS(pDegDiv); k++ )
+      if ( pDivisibleBy2( (pDegDiv->m)[k], resVectorList[j].mon ) )
+        divCount++;
     resVectorList[j].isReduced= (divCount == 1);
+  }
   // create the sets S(k)s
   // a monom x(i)^deg, deg given, is element of the set S(i)
+  // a monom x(i)^deg, deg given, is element of the set S(i)
   // if all x(0)^(polyDegs[0]) ... x(i-1)^(polyDegs[i-1]) DONT divide
   // x(i)^deg and only x(i)^(polyDegs[i]) divides x(i)^deg
   bool doInsert;
+  for ( k= 0; k < iVO->rows(); k++) {
+  for ( k= 0; k < iVO->rows(); k++)
+  {
     //mprPROTInl(" ------------ var:",(*iVO)[k]);
+    for ( j= numVectors - 1; j >= 0; j-- ) {
+    for ( j= numVectors - 1; j >= 0; j-- )
+    {
       //mprPROTPnl("testing monom",resVectorList[j].mon);
+      if ( resVectorList[j].elementOfS == SFREE ) {
+        //mprPROTnl("\tfree");
+        if ( pDivisibleBy2( (pDegDiv->m)[ (*iVO)[k] ], resVectorList[j].mon ) ) {
+          //mprPROTPnl("\tdivisible by ",(pDegDiv->m)[ (*iVO)[k] ]);
+          doInsert=TRUE;
+          for ( i= 0; i < k; i++ ) {
+            //mprPROTPnl("\tchecking db ",(pDegDiv->m)[ (*iVO)[i] ]);
+            if ( pDivisibleBy2( (pDegDiv->m)[ (*iVO)[i] ], resVectorList[j].mon ) ) {
+              //mprPROTPnl("\t and divisible by",(pDegDiv->m)[ (*iVO)[i] ]);
+              doInsert=FALSE;
+              break;
+            }
+          }
+          if ( doInsert ) {
+            //mprPROTInl("\t------------------> S ",(*iVO)[k]);
+            resVectorList[j].elementOfS= (*iVO)[k];
+            resVectorList[j].dividedBy= pCopy( (pDegDiv->m)[ (*iVO)[i] ] );
+            pSetm( resVectorList[j].dividedBy );
+          }
+        }
+      if ( resVectorList[j].elementOfS == SFREE )
+      {
+        //mprPROTnl("\tfree");
+        if ( pDivisibleBy2( (pDegDiv->m)[ (*iVO)[k] ], resVectorList[j].mon ) )
+        {
+          //mprPROTPnl("\tdivisible by ",(pDegDiv->m)[ (*iVO)[k] ]);
+          doInsert=TRUE;
+          for ( i= 0; i < k; i++ )
+          {
+            //mprPROTPnl("\tchecking db ",(pDegDiv->m)[ (*iVO)[i] ]);
+            if ( pDivisibleBy2( (pDegDiv->m)[ (*iVO)[i] ], resVectorList[j].mon ) )
+            {
+              //mprPROTPnl("\t and divisible by",(pDegDiv->m)[ (*iVO)[i] ]);
+              doInsert=FALSE;
+              break;
+            }
+          }
+          if ( doInsert )
+          {
+            //mprPROTInl("\t------------------> S ",(*iVO)[k]);
+            resVectorList[j].elementOfS= (*iVO)[k];
+            resVectorList[j].dividedBy= pCopy( (pDegDiv->m)[ (*iVO)[i] ] );
+            pSetm( resVectorList[j].dividedBy );
+          }
+        }
+      }
+    }
 …
   subSize= 0;
   int sub;
+  for ( i= 0; i < polyDegs->rows(); i++ ) {
+  for ( i= 0; i < polyDegs->rows(); i++ )
+  {
     sub= 1;
     for ( k= 0; k < polyDegs->rows(); k++ )
+    for ( k= 0; k < polyDegs->rows(); k++ )
       if ( i != k ) sub*= (*polyDegs)[k];
     subSize+= sub;
 …
   // Print a list of monoms and their properties
   Print("// \n");
+  for ( j= numVectors - 1; j >= 0; j-- ) {
+  for ( j= numVectors - 1; j >= 0; j-- )
+  {
     Print("// %s, S(%d),  db ",
           resVectorList[j].isReduced?"reduced":"nonreduced",
           resVectorList[j].elementOfS);
     pWrite0(resVectorList[j].dividedBy);
+          resVectorList[j].isReduced?"reduced":"nonreduced",
+          resVectorList[j].elementOfS);
+    pWrite0(resVectorList[j].dividedBy);
     Print("  monom ");
     pWrite(resVectorList[j].mon);
 …
   // make sure that the homogenization variable goes last!
   intvec iVO( pVariables );
+  if ( linPolyS != SNONE ) {
+  if ( linPolyS != SNONE )
+  {
     iVO[pVariables - 1]= linPolyS;
     int p=0;
+    for ( k= pVariables - 1; k >= 0; k-- ) {
+      if ( k != linPolyS ) {
+        iVO[p]= k;
+        p++;
+      }
+    }
+  } else {
+    for ( k= pVariables - 1; k >= 0; k-- )
+    {
+      if ( k != linPolyS )
+      {
+        iVO[p]= k;
+        p++;
+      }
+    }
+  }
+  else
+  {
     linPolyS= 0;
     for ( k= 0; k < pVariables; k++ )
+    for ( k= 0; k < pVariables; k++ )
       iVO[k]= pVariables - k - 1;
+  }
 …
   // generate "matrix"
+  for ( k= numVectors - 1; k >= 0; k-- ) {
+    if ( resVectorList[k].elementOfS != linPolyS ) {
+  for ( k= numVectors - 1; k >= 0; k-- )
+  {
+    if ( resVectorList[k].elementOfS != linPolyS )
+    {
       // column k is a normal column with numerical or symbolic entries
       // init stuff
 …
       // fill in "matrix"
+      while ( pi ) {
+        matEntry= nCopy(pGetCoeff(pi));
+        pmatchPos= pHead0( pi );
+        pSetCoeff( pmatchPos, nInit(1) );
+        pSetm( pmatchPos );
+        for ( i= 0; i < numVectors; i++)
+          if ( pEqual( pmatchPos, resVectorList[i].mon ) )
+            break;
+        resVectorList[k].numColVector[numVectors - i - 1] = nCopy(matEntry);
+        pDelete( &pmatchPos );
+        nDelete( &matEntry );
+        pIter( pi );
+      while ( pi )
+      {
+        matEntry= nCopy(pGetCoeff(pi));
+        pmatchPos= pHead0( pi );
+        pSetCoeff( pmatchPos, nInit(1) );
+        pSetm( pmatchPos );
+        for ( i= 0; i < numVectors; i++)
+          if ( pEqual( pmatchPos, resVectorList[i].mon ) )
+            break;
+        resVectorList[k].numColVector[numVectors - i - 1] = nCopy(matEntry);
+        pDelete( &pmatchPos );
+        nDelete( &matEntry );
+        pIter( pi );
+      }
       pDelete( &pi );
+    } else {
+    }
+    else
+    {
       // column is a special column, i.e. is generated by S0 and F0
       // safe only the positions of the ui's in the column
 …
       j=0;
       while ( pi ) { // fill in "matrix"
         pmp= pMult( pHead( pi ), pCopy( factor ) );
         pSetm( pmp );pTest( pi );
+        for ( i= 0; i < numVectors; i++)
           if ( pEqual( pmp, resVectorList[i].mon ) )
             break;
         resVectorList[k].numColParNr[j]= i;
         pDelete( &pmp );
         pIter( pi );
         j++;
+        pmp= pMult( pHead( pi ), pCopy( factor ) );
+        pSetm( pmp );pTest( pi );
+        for ( i= 0; i < numVectors; i++)
+          if ( pEqual( pmp, resVectorList[i].mon ) )
+            break;
+        resVectorList[k].numColParNr[j]= i;
+        pDelete( &pmp );
+        pIter( pi );
+        j++;
+      }
       pDelete( &pi );
 …
+}
 resVector *resMatrixDense::getMVector(int i)
+{
+resVector *resMatrixDense::getMVector(int i)
+{
   assume( i >= 0 && i < numVectors );
   return &resVectorList[i];
+  return &resVectorList[i];
+}
 const ideal resMatrixDense::getMatrix()
+{
+{
   int k,i,j;
   // copy matrix
   matrix resmat= mpNew(numVectors,numVectors);
+  for (i=1; i <= numVectors; i++) {
+    for (j=1; j <= numVectors; j++ ) {
+      if ( MATELEM(m,i,j) && pGetCoeff(MATELEM(m,i,j)) ) {
+        MATELEM(resmat,i,j)= pCopy( MATELEM(m,i,j) );
+      }
+    }
+  }
+  for (i=0; i < numVectors; i++) {
+    if ( resVectorList[i].elementOfS == linPolyS ) {
+      for (j=1; j <= pVariables; j++ ) {
+        if ( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]) )
+          pDelete( &MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]) );
+        MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1])= pOne();
+        // FIX ME
+        if ( true ) {
+          pSetCoeff( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), nPar(j) );
+        } else {
+          pSetExp( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), j, 1 );
+        }
+  for (i=1; i <= numVectors; i++)
+  {
+    for (j=1; j <= numVectors; j++ )
+    {
+      if ( MATELEM(m,i,j) && pGetCoeff(MATELEM(m,i,j)) )
+      {
+        MATELEM(resmat,i,j)= pCopy( MATELEM(m,i,j) );
+      }
+    }
+  }
+  for (i=0; i < numVectors; i++)
+  {
+    if ( resVectorList[i].elementOfS == linPolyS )
+    {
+      for (j=1; j <= pVariables; j++ )
+      {
+        if ( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]) )
+          pDelete( &MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]) );
+        MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1])= pOne();
+        // FIX ME
+        if ( true )
+        {
+          pSetCoeff( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), nPar(j) );
+        }
+        else
+        {
+          pSetExp( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), j, 1 );
+        }
+      }
+    }
 …
   ideal resmod= idMatrix2Module(resmat);
+  if ( resmat != NULL ) {
+  if ( resmat != NULL )
+  {
     for ( i= 1; i <= numVectors; i++ )
       for ( j= 1; j <= numVectors; j++ )
         pDelete( &MATELEM(resmat , i, j) );
+        pDelete( &MATELEM(resmat , i, j) );
     Free( (ADDRESS)resmat->m, numVectors * numVectors * sizeof(poly) );
     Free( (ADDRESS)resmat, sizeof(ip_smatrix) );
 …
   j=1;
+  for ( k= numVectors - 1; k >= 0; k-- ) {
+  for ( k= numVectors - 1; k >= 0; k-- )
+  {
     vecp= getMVector(k);
     if ( vecp->isReduced ) continue;
     l=1;
+    for ( i= numVectors - 1; i >= 0; i-- ) {
+    for ( i= numVectors - 1; i >= 0; i-- )
+    {
       if ( getMVector(i)->isReduced ) continue;
+      if ( !nIsZero(vecp->getElemNum(numVectors - i - 1)) ) {
+        MATELEM(resmat,j,l)= pInit();
+        MATELEM(resmat,j,l)= pCopy( vecp->getElem(numVectors-i-1) );
+      if ( !nIsZero(vecp->getElemNum(numVectors - i - 1)) )
+      {
+        MATELEM(resmat,j,l)= pInit();
+        MATELEM(resmat,j,l)= pCopy( vecp->getElem(numVectors-i-1) );
+      }
       l++;
+    }
     j++;
+  }
+  }
   ideal resmod= idMatrix2Module(resmat);
+  if ( resmat != NULL ) {
+  if ( resmat != NULL )
+  {
     for ( i= 1; i <= numVectors; i++ )
       for ( j= 1; j <= numVectors; j++ )
         pDelete( &MATELEM(resmat , i, j) );
+        pDelete( &MATELEM(resmat , i, j) );
     Free( (ADDRESS)resmat->m, numVectors * numVectors * sizeof(poly) );
     Free( (ADDRESS)resmat, sizeof(ip_smatrix) );
 …
   // copy evaluation point into matrix
   // p0, p1, ..., pn replace u0, u1, ..., un
+  for ( k= numVectors - 1; k >= 0; k-- ) {
+    if ( linPolyS == getMVector(k)->elementOfS ) {
+      for ( i= 0; i < pVariables; i++ ) {
+        pSetCoeff( MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i]),
+                   nCopy(evpoint[i]) );
+  for ( k= numVectors - 1; k >= 0; k-- )
+  {
+    if ( linPolyS == getMVector(k)->elementOfS )
+    {
+      for ( i= 0; i < pVariables; i++ )
+      {
+        pSetCoeff( MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i]),
+                   nCopy(evpoint[i]) );
+      }
+    }
 …
   // avoid errors for det==0
   number numres;
+  if ( res && pGetCoeff( res ) ) {
+  if ( res && pGetCoeff( res ) )
+  {
     numres= nCopy( pGetCoeff( res ) );
+  } else {
+  }
+  else
+  {
     numres= nInit(0);
     mprPROT("0");
+  }
   pDelete( &res );
+  pDelete( &res );
   mprSTICKYPROT(ST__DET);
 …
   matrix mat= mpNew( subSize, subSize );
+  for ( i= 1; i <= MATROWS( mat ); i++ ) {
+    for ( j= 1; j <= MATCOLS( mat ); j++ ) {
+  for ( i= 1; i <= MATROWS( mat ); i++ )
+  {
+    for ( j= 1; j <= MATCOLS( mat ); j++ )
+    {
       MATELEM(mat,i,j)= pInit();
       pSetCoeff0( MATELEM(mat,i,j), nInit(0) );
 …
+  }
   j=1;
+  for ( k= numVectors - 1; k >= 0; k-- ) {
+  for ( k= numVectors - 1; k >= 0; k-- )
+  {
     vecp= getMVector(k);
     if ( vecp->isReduced ) continue;
     l=1;
+    for ( i= numVectors - 1; i >= 0; i-- ) {
+    for ( i= numVectors - 1; i >= 0; i-- )
+    {
       if ( getMVector(i)->isReduced ) continue;
+      if ( vecp->getElemNum(numVectors - i - 1) && !nIsZero(vecp->getElemNum(numVectors - i - 1)) ) {
+        pSetCoeff(MATELEM(mat, j , l ), nCopy(vecp->getElemNum(numVectors - i - 1)));
+      } /* else {
+           MATELEM(mat, j , l )= pOne();
+           pSetCoeff(MATELEM(mat, j , l ), nInit(0) );
+           }
+        */
+      if ( vecp->getElemNum(numVectors - i - 1) && !nIsZero(vecp->getElemNum(numVectors - i - 1)) )
+      {
+        pSetCoeff(MATELEM(mat, j , l ), nCopy(vecp->getElemNum(numVectors - i - 1)));
+      }
+      /* else
+      {
+           MATELEM(mat, j , l )= pOne();
+           pSetCoeff(MATELEM(mat, j , l ), nInit(0) );
+      }
+      */
       l++;
+    }
     j++;
+  }
+  }
   poly res= singclap_det( mat );
   number numres;
+  if ( res && pGetCoeff( res ) ) {
+  if ( res && pGetCoeff( res ) )
+  {
     numres= nCopy(pGetCoeff( res ));
+    pDelete( &res );
+  } else {
+    pDelete( &res );
+  }
+  else
+  {
     numres= nInit(0);
+  }
 …
   mpz_init(md);mpz_set_ui(md,1);
   mpz_init(mn);mpz_set_ui(mn,1);
   mpz_fac_ui(m,n+d);
   mpz_fac_ui(md,d);
 …
 uResultant::uResultant( const ideal _gls, const resMatType _rmt, BOOLEAN extIdeal )
   : rmt( _rmt )
+{
+  if ( extIdeal ) {
+{
+  if ( extIdeal )
+  {
     // extend given ideal by linear poly F0=u0x0 + u1x1 +...+ unxn
     gls= extendIdeal( _gls, linearPoly( rmt ), rmt );
     n= IDELEMS( gls );
   } else
+  } else
     gls= idCopy( _gls );
+  switch ( rmt ) {
+  switch ( rmt )
+  {
   case sparseResMat:
     resMat= new resMatrixSparse( gls );
 …
     break;
   default:
     Werror("uResultant::uResultant: Unknown resultant matrix type choosen!");
+    WerrorS("uResultant::uResultant: Unknown resultant matrix type choosen!");
+  }
+}
 …
+{
   ideal newGls= idCopy( gls );
   newGls->m= (poly *)ReAlloc( newGls->m,
+                              IDELEMS(gls) * sizeof(poly),
                               (IDELEMS(gls) + 1) * sizeof(poly) );
+  newGls->m= (poly *)ReAlloc( newGls->m,
+                              IDELEMS(gls) * sizeof(poly),
+                              (IDELEMS(gls) + 1) * sizeof(poly) );
   IDELEMS(newGls)++;
+  switch ( rmt ) {
+  switch ( rmt )
+  {
   case sparseResMat:
   case denseResMat:
+    {
       int i;
+      for ( i= IDELEMS(newGls)-1; i > 0; i-- ) {
+        newGls->m[i]= newGls->m[i-1];
+      for ( i= IDELEMS(newGls)-1; i > 0; i-- )
+      {
+        newGls->m[i]= newGls->m[i-1];
+      }
       newGls->m[0]= linPoly;
 …
     break;
   default:
     Werror("uResultant::extendIdeal: Unknown resultant matrix type choosen!");
+  }
+    WerrorS("uResultant::extendIdeal: Unknown resultant matrix type choosen!");
+  }
   return( newGls );
+}
 …
   poly actlp, rootlp= newlp;
+  for ( i= 1; i <= pVariables; i++ ) {
+  for ( i= 1; i <= pVariables; i++ )
+  {
     actlp= newlp;
     pSetExp( actlp, i, 1 );
 …
   pDelete( &newlp );
+  if ( rmt == sparseResMat ) {
+  if ( rmt == sparseResMat )
+  {
     newlp= pOne();
     actlp->next= newlp;
 …
   // long mdg= (facul(tdg+n-1) / facul( tdg )) / facul( n - 1 );
   long mdg= over( n-1, tdg );
   // maxiaml number of terms in a polynom of degree tdg
   long l=(long)pow( tdg+1, n );
 …
 #endif
   // we need mdg results of D(p0,p1,...,pn)
+  // we need mdg results of D(p0,p1,...,pn)
   number *presults;
   presults= (number *)Alloc( mdg * sizeof( number ) );
 …
   // initial evaluatoin point
   p=1;
+  for (i=0; i < n; i++) {
+    // init pevpoint with primes 3,5,7,11, ...
+  for (i=0; i < n; i++)
+  {
+    // init pevpoint with primes 3,5,7,11, ...
     p= nextPrime( p );
     pevpoint[i]=nInit( p );
 …
   // evaluate the determinant in the points pev^0, pev^1, ..., pev^mdg
   mprPROTnl("// evaluating:");
+  for ( i=0; i < mdg; i++ ) {
+    for (j=0; j < n; j++) {
+  for ( i=0; i < mdg; i++ )
+  {
+    for (j=0; j < n; j++)
+    {
       nDelete( &pev[j] );
       nPower(pevpoint[j],i,&pev[j]);
 …
   mprPROTnl("// interpolating:");
   number *ncpoly;
+  {
+  {
     vandermonde vm( mdg, n, tdg, pevpoint );
     ncpoly= vm.interpolateDense( presults );
 …
   if ( subDetVal != NULL ) {   // divide by common factor
       number detdiv;
+      for ( i= 0; i <= mdg; i++ ) {
+        detdiv= nDiv( ncpoly[i], subDetVal );
+        nNormalize( detdiv );
+        nDelete( &ncpoly[i] );
+        ncpoly[i]= detdiv;
+      }
+    }
+      for ( i= 0; i <= mdg; i++ )
+      {
+        detdiv= nDiv( ncpoly[i], subDetVal );
+        nNormalize( detdiv );
+        nDelete( &ncpoly[i] );
+        ncpoly[i]= detdiv;
+      }
+    }
 #ifdef mprDEBUG_ALL
   PrintLn();
+  for ( i=0; i < mdg; i++ ) {
+    nPrint(ncpoly[i]); Print(" --- ");
+  for ( i=0; i < mdg; i++ )
+  {
+    nPrint(ncpoly[i]); Print(" --- ");
+  }
   PrintLn();
 …
   // prepare ncpoly for later use
   number nn=nInit(0);
+  for ( i=0; i < mdg; i++ ) {
+    if ( nEqual(ncpoly[i],nn) ) {
+  for ( i=0; i < mdg; i++ )
+  {
+    if ( nEqual(ncpoly[i],nn) )
+    {
       nDelete( &ncpoly[i] );
       ncpoly[i]=NULL;
 …
   long c=0;
+  for ( i=0; i < l; i++ ) {
+    if ( sum == tdg ) {
+      if ( ncpoly[c] ) {
+        poly p= pOne();
+        if ( rmt == denseResMat ) {
+          for ( j= 0; j < n; j++ ) pSetExp( p, j+1, exp[j] );
+        } else if ( rmt == sparseResMat ) {
+          for ( j= 1; j < n; j++ ) pSetExp( p, j, exp[j] );
+        }
+        pSetCoeff( p, ncpoly[c] );
+        pSetm( p );
+        if (result) result= pAdd( pCopy(result), pCopy(p) );
+        else result= pCopy( p );
+        pDelete( &p );
+  for ( i=0; i < l; i++ )
+  {
+    if ( sum == tdg )
+    {
+      if ( ncpoly[c] )
+      {
+        poly p= pOne();
+        if ( rmt == denseResMat )
+        {
+          for ( j= 0; j < n; j++ ) pSetExp( p, j+1, exp[j] );
+        }
+        else if ( rmt == sparseResMat )
+        {
+          for ( j= 1; j < n; j++ ) pSetExp( p, j, exp[j] );
+        }
+        pSetCoeff( p, ncpoly[c] );
+        pSetm( p );
+        if (result) result= pAdd( pCopy(result), pCopy(p) );
+        else result= pCopy( p );
+        pDelete( &p );
+      }
       c++;
 …
     sum=0;
     exp[0]++;
+    for ( j= 0; j < n - 1; j++ ) {
+      if ( exp[j] > tdg ) {
+        exp[j]= 0;
+        exp[j + 1]++;
+    for ( j= 0; j < n - 1; j++ )
+    {
+      if ( exp[j] > tdg )
+      {
+        exp[j]= 0;
+        exp[j + 1]++;
+      }
       sum+=exp[j];
+    }
+    }
     sum+=exp[n-1];
+  }
+  }
   pSetm( result );
   pTest( result );
   return result;
+}
 …
   // evaluate D in tdg+1 distinct points, so
   // we need tdg+1 results of D(p0,1,0,...,0) =
+  // we need tdg+1 results of D(p0,1,0,...,0) =
   //              c(0)*u0^tdg + c(1)*u0^tdg-1 + ... + c(tdg-1)*u0 + c(tdg)
   number *presults;
 …
   // this gives us n-1 evaluations
   p=3;
+  for ( uvar= 0; uvar < loops; uvar++ ) {
+  for ( uvar= 0; uvar < loops; uvar++ )
+  {
     // generate initial evaluation point
+    if ( matchUp ) {
+      for (i=0; i < n; i++) {
+        // prime(random number) between 1 and MAXEVPOINT
+        nDelete( &pevpoint[i] );
+        if ( i == 0 ) {
+          //p= nextPrime( p );
+          pevpoint[i]= nInit( p );
+        } else
+          if ( i <= uvar + 2 ) {
+            pevpoint[i]=nInit(IsPrime(1+(int) (MAXEVPOINT*rand()/(RAND_MAX+1.0))));
+            //pevpoint[i]=nInit(383);
+          } else pevpoint[i]=nInit(0);
+        mprPROTNnl(" ",pevpoint[i]);
+      }
+    } else {
+      for (i=0; i < n; i++) {
+        // init pevpoint with  prime,0,...0,1,0,...,0
+        nDelete( &pevpoint[i] );
+        if ( i == 0 ) {
+          //p=nextPrime( p );
+          pevpoint[i]=nInit( p );
+        }
+        else
+          if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
+          else pevpoint[i]= nInit(0);
+        mprPROTNnl(" ",pevpoint[i]);
+    if ( matchUp )
+    {
+      for (i=0; i < n; i++)
+      {
+        // prime(random number) between 1 and MAXEVPOINT
+        nDelete( &pevpoint[i] );
+        if ( i == 0 )
+        {
+          //p= nextPrime( p );
+          pevpoint[i]= nInit( p );
+        }
+        else if ( i <= uvar + 2 )
+        {
+          pevpoint[i]=nInit(IsPrime(1+(int) (MAXEVPOINT*rand()/(RAND_MAX+1.0))));
+          //pevpoint[i]=nInit(383);
+        }
+        else
+          pevpoint[i]=nInit(0);
+        mprPROTNnl(" ",pevpoint[i]);
+      }
+    }
+    else
+    {
+      for (i=0; i < n; i++)
+      {
+        // init pevpoint with  prime,0,...0,1,0,...,0
+        nDelete( &pevpoint[i] );
+        if ( i == 0 )
+        {
+          //p=nextPrime( p );
+          pevpoint[i]=nInit( p );
+        }
+        else
+          if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
+          else pevpoint[i]= nInit(0);
+        mprPROTNnl(" ",pevpoint[i]);
+      }
+    }
     // prepare aktual evaluation point
+    for (i=0; i < n; i++) {
+    for (i=0; i < n; i++)
+    {
       nDelete( &pev[i] );
       pev[i]= nCopy( pevpoint[i] );
+    }
     // evaluate the determinant in the points pev^0, pev^1, ..., pev^tdg
+    for ( i=0; i <= tdg; i++ ) {
+    for ( i=0; i <= tdg; i++ )
+    {
       nDelete( &pev[0] );
       nPower(pevpoint[0],i,&pev[0]);          // new evpoint
 …
     mprSTICKYPROT("\n");
     // now interpolate
+    // now interpolate
     vandermonde vm( tdg + 1, 1, tdg, pevpoint, FALSE );
     number *ncpoly= vm.interpolateDense( presults );
 …
     if ( subDetVal != NULL ) {  // divide by common factor
       number detdiv;
+      for ( i= 0; i <= tdg; i++ ) {
+        detdiv= nDiv( ncpoly[i], subDetVal );
+        nNormalize( detdiv );
+        nDelete( &ncpoly[i] );
+        ncpoly[i]= detdiv;
+      for ( i= 0; i <= tdg; i++ )
+      {
+        detdiv= nDiv( ncpoly[i], subDetVal );
+        nNormalize( detdiv );
+        nDelete( &ncpoly[i] );
+        ncpoly[i]= detdiv;
+      }
+    }
 …
 #ifdef mprDEBUG_ALL
     PrintLn();
+    for ( i=0; i <= tdg; i++ ) {
+      nPrint(ncpoly[i]); Print(" --- ");
+    for ( i=0; i <= tdg; i++ )
+    {
+      nPrint(ncpoly[i]); Print(" --- ");
+    }
     PrintLn();
 …
     // save results
     roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
                                 (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
                                 loops );
+    roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
+                                (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
+                                loops );
+  }
 …
   // or D(u0,k1,k2,0,...,0), D(u0,k1,k2,k3,0,...,0), ..., D(u0,k1,k2,k3,...,kn)
   p=3;
+  for ( uvar= 0; uvar < loops; uvar++ ) {
+  for ( uvar= 0; uvar < loops; uvar++ )
+  {
     // generate initial evaluation point
+    if ( matchUp ) {
+      for (i=0; i < n; i++) {
+        // prime(random number) between 1 and MAXEVPOINT
+        nDelete( &pevpoint[i] );
+        if ( i <= uvar + 2 ) {
+          pevpoint[i]=nInit(IsPrime(1+(int) (MAXEVPOINT*rand()/(RAND_MAX+1.0))));
+          //pevpoint[i]=nInit(383);
+        } else pevpoint[i]=nInit(0);
+        mprPROTNnl(" ",pevpoint[i]);
+      }
+    } else {
+      for (i=0; i < n; i++) {
+        // init pevpoint with  prime,0,...0,-1,0,...,0
+        nDelete( &(pevpoint[i]) );
+        if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
+        else pevpoint[i]= nInit(0);
+        mprPROTNnl(" ",pevpoint[i]);
+    if ( matchUp )
+    {
+      for (i=0; i < n; i++)
+      {
+        // prime(random number) between 1 and MAXEVPOINT
+        nDelete( &pevpoint[i] );
+        if ( i <= uvar + 2 )
+        {
+          pevpoint[i]=nInit(IsPrime(1+(int) (MAXEVPOINT*rand()/(RAND_MAX+1.0))));
+          //pevpoint[i]=nInit(383);
+        } else pevpoint[i]=nInit(0);
+        mprPROTNnl(" ",pevpoint[i]);
+      }
+    }
+    else
+    {
+      for (i=0; i < n; i++)
+      {
+        // init pevpoint with  prime,0,...0,-1,0,...,0
+        nDelete( &(pevpoint[i]) );
+        if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
+        else pevpoint[i]= nInit(0);
+        mprPROTNnl(" ",pevpoint[i]);
+      }
+    }
 …
     piter= pures;
+    for ( i= tdg; i >= 0; i-- ) {
+    for ( i= tdg; i >= 0; i-- )
+    {
       //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
+      if ( piter && pTotaldegree(piter) == i ) {
+        ncpoly[i]= nCopy( pGetCoeff( piter ) );
+        pIter( piter );
+      } else {
+        ncpoly[i]= nInit(0);
+      if ( piter && pTotaldegree(piter) == i )
+      {
+        ncpoly[i]= nCopy( pGetCoeff( piter ) );
+        pIter( piter );
+      }
+      else
+      {
+        ncpoly[i]= nInit(0);
+      }
       mprPROTNnl("", ncpoly[i] );
+    }
+    }
     mprSTICKYPROT(ST_BASE_EV); // .
     if ( subDetVal != NULL ) {  // divide by common factor
       number detdiv;
+      for ( i= 0; i <= tdg; i++ ) {
+        detdiv= nDiv( ncpoly[i], subDetVal );
+        nNormalize( detdiv );
+        nDelete( &ncpoly[i] );
+        ncpoly[i]= detdiv;
+      for ( i= 0; i <= tdg; i++ )
+      {
+        detdiv= nDiv( ncpoly[i], subDetVal );
+        nNormalize( detdiv );
+        nDelete( &ncpoly[i] );
+        ncpoly[i]= detdiv;
+      }
+    }
 …
     // save results
     roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
                                 (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
                                 loops );
+    roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
+                                (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
+                                loops );
+  }
 …
   i+=2;
   int j= IsPrime( i );
+  while ( j <= init ) {
+  while ( j <= init )
+  {
     i+=2;
     j= IsPrime( i );
 …
 // compile-command-1: "make installg" ***
 // compile-command-2: "make install" ***
 // End: ***
+// End: ***
 // in folding: C-c x

Singular/mpr_base.h

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: mpr_base.h,v 1.1 1999-06-28 12:48:12 wenk Exp $ */
+/* $Id: mpr_base.h,v 1.2 1999-06-28 16:06:25 Singular Exp $ */
 /*
+/*
 * ABSTRACT - multipolynomial resultants - resultant matrices
 *            ( sparse, dense, u-resultant solver )
 …
  * Base class for sparse and dense u-Resultant computation
  */
+class resMatrixBase {
+class resMatrixBase
+{
 public:
   /* state of the resultant */
 …
   virtual const number getSubDet() { return NULL; }
   virtual const long getDetDeg() const { return totDeg; }
+  virtual const long getDetDeg() const { return totDeg; }
   virtual const IStateType initState() const { return istate; }
+  virtual const IStateType initState() const { return istate; }
 protected:
 …
   int totDeg;
 private:
   /* disables the copy constructor */
 …
  * Base class for solving 0-dim poly systems using u-resultant
  */
+class uResultant {
+class uResultant
+{
 public:
   enum resMatType { none, sparseResMat, denseResMat };
 …
   rootContainer ** specializeInU( BOOLEAN matchUp= false, const number subDetVal= NULL );
   resMatrixBase * accessResMat() { return resMat; }
+  resMatrixBase * accessResMat() { return resMat; }
 private:
   /* deactivated copy constructor */
   uResultant( const uResultant & );
   ideal extendIdeal( const ideal gls, poly linPoly, const resMatType rmt );
   poly linearPoly( const resMatType rmt );
 …
 // compile-command-2: "make install" ***
 // compile-command: "make installg" ***
+// End: ***
+// End: ***

Singular/mpr_complex.cc

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: mpr_complex.cc,v 1.7 1999-06-28 12:48:12 wenk Exp $ */
+/* $Id: mpr_complex.cc,v 1.8 1999-06-28 16:06:26 Singular Exp $ */
 /*
 …
   gmp_float r;
+  if ( rField_is_Q() ) {
+  if ( rField_is_Q() )
+  {
     if ( num != NULL )
+    {
+      if (SR_HDL(num) & SR_INT)
+      {
+        if (SR_HDL(num) & SR_INT)
+          {
+            r= SR_TO_INT(num);
+          }
+        else
+          {
+            if ( num->s == 0 )
+              {
+                nlNormalize( num );
+              }
+            if (SR_HDL(num) & SR_INT)
+              {
+                r= SR_TO_INT(num);
+              }
+            else
+              {
+                if ( num->s != 3 )
+                  {
+                    r= &num->z;
+                    r/= (gmp_float)&num->n;
+                  }
+                else
+                  {
+                    r= &num->z;
+                  }
+              }
+          }
+        r= SR_TO_INT(num);
+      }
+      else
+      {
+        if ( num->s == 0 )
+        {
+          nlNormalize( num );
+        }
+        if (SR_HDL(num) & SR_INT)
+        {
+          r= SR_TO_INT(num);
+        }
+        else
+        {
+          if ( num->s != 3 )
+          {
+            r= &num->z;
+            r/= (gmp_float)&num->n;
+          }
+          else
+          {
+            r= &num->z;
+          }
+        }
+      }
+    }
     else
+      {
+        r= 0.0;
+      }
+  } else if (rField_is_long_R() || rField_is_long_C()) {
+    {
+      r= 0.0;
+    }
+  }
+  else if (rField_is_long_R() || rField_is_long_C())
+  {
     r= *(gmp_float*)num;
+  } else if ( rField_is_R() ) {
+  }
+  else if ( rField_is_R() )
+  {
     // Add some code here :-)
+    Werror("Ground field not implemented!");
+  } else {
+    Werror("Ground field not implemented!");
+    WerrorS("Ground field not implemented!");
+  }
+  else
+  {
+    WerrorS("Ground field not implemented!");
+  }
 …
     in_real=floatToStr( c.real(), oprec );         // get real part
     in_imag=floatToStr( abs(c.imag()), oprec );    // get imaginary part
+    if (rField_is_long_C()) {
+    if (rField_is_long_C())
+    {
       out=(char*)AllocL((strlen(in_real)+strlen(in_imag)+5+strlen(currRing->parameter[0]))*sizeof(char));
       sprintf(out,"%s%s%s*%s",in_real,c.imag().sign()>=0?"+":"-",currRing->parameter[0],in_imag);
+    } else {
+    }
+    else
+    {
       out=(char*)AllocL( (strlen(in_real)+strlen(in_imag)+8) * sizeof(char));
       sprintf(out,"%s%s%s",in_real,c.imag().sign()>=0?" + I ":" - I ",in_imag);

Singular/mpr_complex.h

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: mpr_complex.h,v 1.5 1999-06-28 12:48:13 wenk Exp $ */
 /*
+/* $Id: mpr_complex.h,v 1.6 1999-06-28 16:06:26 Singular Exp $ */
+/*
 * ABSTRACT - multipolynomial resultants - real floating-point numbers using gmp
 *            and complex numbers based on pairs of real floating-point numbers
+*
+*
 */
 //-> include & define stuff
 // must have gmp version >= 2
 extern "C" {
+extern "C" {
 #include <gmp.h>
+}
 …
   inline bool setFromStr( char * in );
   // access
+  // access
   inline const mpf_t *mpfp() const;
 …
 public:
   static void setPrecision( const unsigned long int prec ) {
     gmp_default_prec_bits= prec;
+  }
   static void setEqualBits( const unsigned long int prec ) {
     gmp_needequal_bits= prec;
+  }
   static const unsigned long int getPrecision() {
     return gmp_default_prec_bits;
+  }
   static const unsigned long int getEqualBits() {
     return gmp_needequal_bits;
+  static void setPrecision( const unsigned long int prec ) {
+    gmp_default_prec_bits= prec;
+  }
+  static void setEqualBits( const unsigned long int prec ) {
+    gmp_needequal_bits= prec;
+  }
+  static const unsigned long int getPrecision() {
+    return gmp_default_prec_bits;
+  }
+  static const unsigned long int getEqualBits() {
+    return gmp_needequal_bits;
+  }
 …
 //-> class gmp_complex
 /**
  * @short gmp_complex numbers based on
+ * @short gmp_complex numbers based on
  */
 class gmp_complex
 …
   // gmp_complex <operator> real
   inline friend gmp_complex operator + ( const gmp_complex & a, const gmp_float b_d );
+  inline friend gmp_complex operator + ( const gmp_complex & a, const gmp_float b_d );
   inline friend gmp_complex operator - ( const gmp_complex & a, const gmp_float b_d );
   inline friend gmp_complex operator * ( const gmp_complex & a, const gmp_float b_d );
 …
 // compile-command-1: "make installg" ***
 // compile-command-2: "make install" ***
+// End: ***
+// End: ***

Singular/mpr_global.h

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: mpr_global.h,v 1.3 1999-06-28 12:48:14 wenk Exp $ */
+/* $Id: mpr_global.h,v 1.4 1999-06-28 16:06:27 Singular Exp $ */
 /*
 * ABSTRACT - multipolynomial resultants -
+/*
+* ABSTRACT - multipolynomial resultants -
 *                                global definitions and debugging stuff
 */
 …
 #define mprPROT(msg) Print("%s",msg)
 #define mprPROTnl(msg) Print("%s\n",msg)
 #define mprPROTP(msg,poly) Print("%s",msg);pWrite0(poly)
 #define mprPROTPnl(msg,poly) Print("%s",msg);pWrite(poly)
+#define mprPROTP(msg,poly) Print("%s",msg);pWrite0(poly)
+#define mprPROTPnl(msg,poly) Print("%s",msg);pWrite(poly)
 #define mprPROTI(msg,intval) Print("%s%d",msg,intval)
 #define mprPROTInl(msg,intval) Print("%s%d\n",msg,intval)
 #define mprPROTN(msg,nval) Print("%s",msg);nPrint(nval);
 #define mprPROTNnl(msg,nval) Print("%s",msg);nPrint(nval);PrintLn();
 #else
 #define mprPROT(msg)
+#else
+#define mprPROT(msg)
 #define mprPROTnl(msg)
 #define mprPROTP(msg,poly)
 …
 // compile-command-1: "make installg" ***
 // compile-command-2: "make install" ***
+// End: ***
+// End: ***

Singular/mpr_inout.cc

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: mpr_inout.cc,v 1.1 1999-06-28 12:48:14 wenk Exp $ */
 /*
+/* $Id: mpr_inout.cc,v 1.2 1999-06-28 16:06:28 Singular Exp $ */
+/*
 * ABSTRACT - multipolynomial resultant
 */
 …
 #include "tok.h"
 #include "structs.h"
 #include "subexpr.h"
+#include "subexpr.h"
 #include "polys.h"
 #include "ideals.h"
 …
 #define TIMING_EPR(t,msg) TIMING_END_AND_PRINT(t,msg);TIMING_RESET(t);
+enum mprState{
+    mprOk,
+enum mprState
+{
+    mprOk,
     mprWrongRType,
     mprHasOne,
+    mprHasOne,
     mprInfNumOfVars,
     mprNotReduced,
 …
 //<-
 //-> nPrint(number n)
+//-> nPrint(number n)
 void nPrint(number n)
+{
 …
 void mprPrintError( mprState state, const char * name )
+{
+  switch (state) {
+  switch (state)
+  {
   case mprWrongRType:
     Werror("Unknown resultant matrix type choosen!");
+    WerrorS("Unknown resultant matrix type choosen!");
     break;
   case mprHasOne:
 …
   case mprInfNumOfVars:
     Werror("Numer of elements in given ideal %s must be equal to %d!",
            name,pVariables+1);
+           name,pVariables+1);
     break;
   case mprNotZeroDim:
 …
   case mprNotHomog:
     Werror("The given ideal %s has to be homogeneous in"
            " the first ring variable!",name);
+           " the first ring variable!",name);
     break;
   case mprNotReduced:
 …
     break;
   case mprUnSupField:
     Werror("Ground field not implemented!");
+    WerrorS("Ground field not implemented!");
     break;
   default:
 …
 //-> mprState mprIdealCheck()
 mprState mprIdealCheck( const ideal theIdeal,
+                        const char * name,
+                        uResultant::resMatType mtype,
                         BOOLEAN rmatrix= false )
+mprState mprIdealCheck( const ideal theIdeal,
+                        const char * name,
+                        uResultant::resMatType mtype,
+                        BOOLEAN rmatrix= false )
+{
   mprState state = mprOk;
   int power;
   int k;
+  int k;
   int numOfVars= mtype == uResultant::denseResMat?pVariables-1:pVariables;
   if ( rmatrix ) numOfVars++;
   if ( mtype == uResultant::none )
+  if ( mtype == uResultant::none )
     state= mprWrongRType;
 …
     state= mprInfNumOfVars;
+  for ( k= IDELEMS(theIdeal) - 1; (state == mprOk) && (k >= 0); k-- ) {
+  for ( k= IDELEMS(theIdeal) - 1; (state == mprOk) && (k >= 0); k-- )
+  {
     poly p = (theIdeal->m)[k];
     if ( pIsConstant(p) ) state= mprHasOne;
     else
     if ( (mtype == uResultant::denseResMat) && !pIsHomogeneous(p) )
+    else
+    if ( (mtype == uResultant::denseResMat) && !pIsHomogeneous(p) )
       state=mprNotHomog;
+  }
   if ( !(rField_is_R()||
          rField_is_Q()||
          rField_is_long_R()||
          rField_is_long_C()||
+         (rmatrix && rField_is_Q_a())) )
+         rField_is_Q()||
+         rField_is_long_R()||
+         rField_is_long_C()||
+         (rmatrix && rField_is_Q_a())) )
     state= mprUnSupField;
 …
 uResultant::resMatType determineMType( int imtype )
+{
+  switch ( imtype ) {
+  switch ( imtype )
+  {
   case MPR_DENSE:
     return uResultant::denseResMat;
 …
   res->rtyp = LIST_CMD;
   res->data= (void *)emptylist;
   TIMING_START(mpr_overall);
   // check input ideal ( = polynomial system )
+  if ( mprIdealCheck( gls, arg1->Name(), mtype ) != mprOk ) {
+  if ( mprIdealCheck( gls, arg1->Name(), mtype ) != mprOk )
+  {
     return TRUE;
+  }
 …
   TIMING_START(mpr_constr);
   ures= new uResultant( gls, mtype );
+  if ( ures->accessResMat()->initState() != resMatrixBase::ready ) {
+    Werror("Error occurred during matrix setup!");
+  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
+  {
+    WerrorS("Error occurred during matrix setup!");
     return TRUE;
+  }
 …
   // if dense resultant, check if minor nonsingular
   TIMING_START(mpr_check);
+  if ( mtype == uResultant::denseResMat ) {
+  if ( mtype == uResultant::denseResMat )
+  {
     smv= ures->accessResMat()->getSubDet();
 #ifdef mprDEBUG_PROT
     Print("// Determinant of submatrix: ");nPrint(smv);PrintLn();
 #endif
+    if ( nIsZero(smv) ) {
+      Werror("Unsuitable input ideal: Minor of resultant matrix is singular!");
+    if ( nIsZero(smv) )
+    {
+      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
       return TRUE;
+    }
 …
   TIMING_START(mpr_ures);
   if ( interpolate_det )
     iproots= ures->interpolateDenseSP( false, smv );
   else
+    iproots= ures->interpolateDenseSP( false, smv );
+  else
     iproots= ures->specializeInU( false, smv );
   TIMING_EPR(mpr_ures,   "interpolation ures\t")
 …
   // main task 3: Interpolate specialized resultant polynomials
   TIMING_START(mpr_mures);
   if ( interpolate_det )
+  if ( interpolate_det )
     muiproots= ures->interpolateDenseSP( true, smv );
   else
+  else
     muiproots= ures->specializeInU( true, smv );
   TIMING_EPR(mpr_mures,  "interpolation mures\t")
 …
   arranger->solve_all();
   TIMING_EPR(mpr_solver, "solver time\t\t");
   // get list of roots
+  if ( arranger->success() ) {
+  if ( arranger->success() )
+  {
     TIMING_START(mpr_arrange);
     arranger->arrange();
     TIMING_EPR(mpr_arrange, "arrange time\t\t");
     listofroots= arranger->listOfRoots( gmp_output_digits );
+  } else {
+    Werror("Solver was unable to find any root!");
+  }
+  else
+  {
+    WerrorS("Solver was unable to find any root!");
     return TRUE;
+  }
 …
   emptylist->Clean();
   //Free( (ADDRESS) emptylist, sizeof(slists) );
   TIMING_EPR(mpr_overall,"overall time\t\t")
 …
   // check input ideal ( = polynomial system )
+  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk ) {
+  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
+  {
     return TRUE;
+  }
 …
   delete resMat;
   return FALSE;
+}
 …
   gls= (poly)(arg1->Data());
   int howclean= (int)arg2->Data();
   int deg= pTotaldegree( gls );
   //  int deg= pDeg( gls );
 …
   if ( !(rField_is_R() ||
+         rField_is_Q() ||
+         rField_is_long_R() ||
+         rField_is_long_C()) ) {
+    Werror("Ground field not implemented!\n");
+    return TRUE;
+  }
+  if ( pVariables > 1 ) {
+         rField_is_Q() ||
+         rField_is_long_R() ||
+         rField_is_long_C()) )
+         {
+    WerrorS("Ground field not implemented!");
+    return TRUE;
+  }
+  if ( pVariables > 1 )
+  {
     piter= gls;
+    for ( i= 1; i <= pVariables; i++ )
+      if ( pGetExp( piter, i ) ) {
+        vpos= i;
+        break;
+    for ( i= 1; i <= pVariables; i++ )
+      if ( pGetExp( piter, i ) )
+      {
+        vpos= i;
+        break;
+      }
+    while ( piter ) {
+      for ( i= 1; i <= pVariables; i++ )
+        if ( (vpos != i) && (pGetExp( piter, i ) != 0) ) {
+          Werror("The polynomial %s must be univariate!",arg1->Name());
+          return TRUE;
+        }
+    while ( piter )
+    {
+      for ( i= 1; i <= pVariables; i++ )
+        if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
+        {
+          Werror("The polynomial %s must be univariate!",arg1->Name());
+          return TRUE;
+        }
       pIter( piter );
+    }
 …
   number * pcoeffs= (number *)Alloc( (deg+1) * sizeof( number ) );
   piter= gls;
+  for ( i= deg; i >= 0; i-- ) {
+  for ( i= deg; i >= 0; i-- )
+  {
     //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
+    if ( piter && pTotaldegree(piter) == i ) {
+    if ( piter && pTotaldegree(piter) == i )
+    {
       pcoeffs[i]= nCopy( pGetCoeff( piter ) );
       //nPrint( pcoeffs[i] );Print("  ");
       pIter( piter );
+    } else {
+    }
+    else
+    {
       pcoeffs[i]= nInit(0);
+    }
+  }
+  }
 #ifdef mprDEBUG_PROT
+  for (i=deg; i >= 0; i--) {
+  for (i=deg; i >= 0; i--)
+  {
     nPrint( pcoeffs[i] );Print("  ");
+  }
 …
   rlist->Init( elem );
+  if (rField_is_long_C()) {
+    for ( j= 0; j < elem; j++ ) {
+  if (rField_is_long_C())
+  {
+    for ( j= 0; j < elem; j++ )
+    {
       rlist->m[j].rtyp=NUMBER_CMD;
       rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
       //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
+    }
+  } else {
+    for ( j= 0; j < elem; j++ ) {
+  }
+  else
+  {
+    for ( j= 0; j < elem; j++ )
+    {
       dummy = complexToStr( (*roots)[j], gmp_output_digits );
       rlist->m[j].rtyp=STRING_CMD;
 …
   // p can be a vector of numbers (multivariate polynom)
   //   or one number (univariate polynom)
   // tdg = deg(f)
   int n= IDELEMS( p );
+  // tdg = deg(f)
+  int n= IDELEMS( p );
   int m= IDELEMS( w );
   int tdg= (int)arg3->Data();
   res->data= (void*)NULL;
   // check the input
+  if ( tdg < 1 ) {
+    Werror("Last input parameter must be > 0!");
+    return TRUE;
+  }
+  if ( n != pVariables ) {
+    Werror("Size of first input ideal must be equal to %d!\n",pVariables);
+    return TRUE;
+  }
+  if ( m != (int)pow((double)tdg+1,(int)n) ) {
+    Werror("Size of second input ideal must be equal to %d!\n",(int)pow((double)tdg+1,(int)n));
+  if ( tdg < 1 )
+  {
+    WerrorS("Last input parameter must be > 0!");
+    return TRUE;
+  }
+  if ( n != pVariables )
+  {
+    Werror("Size of first input ideal must be equal to %d!",pVariables);
+    return TRUE;
+  }
+  if ( m != (int)pow((double)tdg+1,(int)n) )
+  {
+    Werror("Size of second input ideal must be equal to %d!",
+      (int)pow((double)tdg+1,(int)n));
     return TRUE;
+  }
   if ( !(rField_is_Q() /* ||
+         rField_is_R() || rField_is_long_R() ||
+         rField_is_long_C()*/ ) ) {
+    Werror("Ground field not implemented!\n");
+         rField_is_R() || rField_is_long_R() ||
+         rField_is_long_C()*/ ) )
+         {
+    WerrorS("Ground field not implemented!");
     return TRUE;
+  }
 …
   number tmp;
   number *pevpoint= (number *)Alloc( n * sizeof( number ) );
+  for ( i= 0; i < n; i++ ) {
+  for ( i= 0; i < n; i++ )
+  {
     pevpoint[i]=nInit(0);
+    if (  (p->m)[i] ) {
+    if (  (p->m)[i] )
+    {
       tmp = pGetCoeff( (p->m)[i] );
+      if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) ) {
+        Free( (ADDRESS)pevpoint, n * sizeof( number ) );
+        Werror("Elements of first input ideal must not be equal to -1, 0, 1!");
+        return TRUE;
+      if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
+      {
+        Free( (ADDRESS)pevpoint, n * sizeof( number ) );
+        WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
+        return TRUE;
+      }
     } else tmp= NULL;
+    if ( !nIsZero(tmp) ) {
+      if ( !pIsConstant((p->m)[i])) {
+        Free( (ADDRESS)pevpoint, n * sizeof( number ) );
+        Werror("Elements of first input ideal must be numbers!");
+        return TRUE;
+    if ( !nIsZero(tmp) )
+    {
+      if ( !pIsConstant((p->m)[i]))
+      {
+        Free( (ADDRESS)pevpoint, n * sizeof( number ) );
+        WerrorS("Elements of first input ideal must be numbers!");
+        return TRUE;
+      }
       pevpoint[i]= nCopy( tmp );
 …
   number *wresults= (number *)Alloc( m * sizeof( number ) );
+  for ( i= 0; i < m; i++ ) {
+  for ( i= 0; i < m; i++ )
+  {
     wresults[i]= nInit(0);
+    if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) ) {
+      if ( !pIsConstant((w->m)[i])) {
+        Free( (ADDRESS)pevpoint, n * sizeof( number ) );
+        Free( (ADDRESS)wresults, m * sizeof( number ) );
+        Werror("Elements of second input ideal must be numbers!");
+        return TRUE;
+    if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
+    {
+      if ( !pIsConstant((w->m)[i]))
+      {
+        Free( (ADDRESS)pevpoint, n * sizeof( number ) );
+        Free( (ADDRESS)wresults, m * sizeof( number ) );
+        WerrorS("Elements of second input ideal must be numbers!");
+        return TRUE;
+      }
       wresults[i]= nCopy(pGetCoeff((w->m)[i]));
 …
 //<-
 //-> function u_resultant_det
+//-> function u_resultant_det
 poly u_resultant_det( ideal gls, int imtype )
+{
 …
   // check input ideal ( = polynomial system )
+  if ( mprIdealCheck( gls, "", mtype ) != mprOk ) {
+  if ( mprIdealCheck( gls, "", mtype ) != mprOk )
+  {
     return emptypoly;
+  }
 …
   // if dense resultant, check if minor nonsingular
+  if ( mtype == uResultant::denseResMat ) {
+  if ( mtype == uResultant::denseResMat )
+  {
     smv= ures->accessResMat()->getSubDet();
 #ifdef mprDEBUG_PROT
     Print("// Determinant of submatrix: ");nPrint(smv); PrintLn();
 #endif
+    if ( nIsZero(smv) ) {
+      Werror("Unsuitable input ideal: Minor of resultant matrix is singular!");
+    if ( nIsZero(smv) )
+    {
+      WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
       return emptypoly;
+    }
 …
 // compile-command-1: "make installg" ***
 // compile-command-2: "make install" ***
 // End: ***
+// End: ***
 // in folding: C-c x

Singular/mpr_inout.h

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: mpr_inout.h,v 1.1 1999-06-28 12:48:15 wenk Exp $ */
+/* $Id: mpr_inout.h,v 1.2 1999-06-28 16:06:28 Singular Exp $ */
 /*
+/*
 * ABSTRACT - multipolynomial resultants - interface to Singular
+*
+*
 */
 …
 /** solve a multipolynomial system using the u-resultant
  * Input ideal must be 0-dimensional and pVariables == IDELEMS(ideal).
  * Resultant method can be MPR_DENSE, which uses Macaulay Resultant (good for
  * dense homogeneous polynoms) or MPR_SPARSE, which uses Sparse Resultant
+ * Resultant method can be MPR_DENSE, which uses Macaulay Resultant (good for
+ * dense homogeneous polynoms) or MPR_SPARSE, which uses Sparse Resultant
  * (Gelfand, Kapranov, Zelevinsky).
  * If interpolate == true then the determinant of the u-resultant will be
  * numerically interpolatet using a Vandermonde System.
+ * If interpolate == true then the determinant of the u-resultant will be
+ * numerically interpolatet using a Vandermonde System.
  * Otherwise, the Sparse Bareiss will be used (faster!).
  * Returns a list containing the roots of the system.
 …
 /** find the (complex) roots an univariate polynomial
  * Determines the roots of an univariate polynomial using Laguerres'
+ * Determines the roots of an univariate polynomial using Laguerres'
  * root-solver. Good for polynomials with low and middle degree (<40).
  * Returns a list containing the roots of the polynomial.
 …
 // compile-command-1: "make installg" ***
 // compile-command-2: "make install" ***
 // End: ***
+// End: ***

Singular/mpr_numeric.cc

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: mpr_numeric.cc,v 1.1 1999-06-28 12:48:16 wenk Exp $ */
 /*
+/* $Id: mpr_numeric.cc,v 1.2 1999-06-28 16:06:29 Singular Exp $ */
+/*
 * ABSTRACT - multipolynomial resultants - numeric stuff
 *            ( root finder, vandermonde system solver, simplex )
+*            ( root finder, vandermonde system solver, simplex )
 */
 …
 //-> vandermonde::*
 vandermonde::vandermonde( const long _cn, const long _n, const long _maxdeg,
                           number *_p, const bool _homog )
+vandermonde::vandermonde( const long _cn, const long _n, const long _maxdeg,
+                          number *_p, const bool _homog )
   : n(_n), cn(_cn), maxdeg(_maxdeg), p(_p), homog(_homog)
+{
 …
   for ( j= 0; j < n; j++ ) exp[j]=0;
+  for ( i= 0; i < l; i++ ) {
+    if ( !homog || (sum == maxdeg) ) {
+      for ( j= 0; j < n; j++ ) {
+        nPower( p[j], exp[j], &tmp );
+        tmp1 = nMult( tmp, x[c] );
+        x[c]= tmp1;
+        nDelete( &tmp );
+  for ( i= 0; i < l; i++ )
+  {
+    if ( !homog || (sum == maxdeg) )
+    {
+      for ( j= 0; j < n; j++ )
+      {
+        nPower( p[j], exp[j], &tmp );
+        tmp1 = nMult( tmp, x[c] );
+        x[c]= tmp1;
+        nDelete( &tmp );
+      }
       c++;
 …
     exp[0]++;
     sum=0;
+    for ( j= 0; j < n - 1; j++ ) {
+      if ( exp[j] > maxdeg ) {
+        exp[j]= 0;
+        exp[j + 1]++;
+        }
+    for ( j= 0; j < n - 1; j++ )
+    {
+      if ( exp[j] > maxdeg )
+      {
+        exp[j]= 0;
+        exp[j + 1]++;
+        }
       sum+= exp[j];
+    }
 …
   for ( j= 0; j < n+1; j++ ) exp[j]=0;
+  for ( i= 0; i < l; i++ ) {
+    if ( (!homog || (sum == maxdeg)) && q[i] && !nIsZero(q[i]) ) {
+  for ( i= 0; i < l; i++ )
+  {
+    if ( (!homog || (sum == maxdeg)) && q[i] && !nIsZero(q[i]) )
+    {
       pnew= pOne();
       pSetCoeff(pnew,q[i]);
       pSetExpV(pnew,exp);
+      if ( pit ) {
+        pNext(pnew)= pit;
+        pit= pnew;
+      } else {
+        pit= pnew;
+        pNext(pnew)= NULL;
+      }
+      if ( pit )
+      {
+        pNext(pnew)= pit;
+        pit= pnew;
+      }
+      else
+      {
+        pit= pnew;
+        pNext(pnew)= NULL;
+      }
       pSetm(pit);
+    }
     exp[1]++;
     sum=0;
+    for ( j= 1; j < n; j++ ) {
+      if ( exp[j] > maxdeg ) {
+        exp[j]= 0;
+        exp[j + 1]++;
+        }
+    for ( j= 1; j < n; j++ )
+    {
+      if ( exp[j] > maxdeg )
+      {
+        exp[j]= 0;
+        exp[j + 1]++;
+        }
       sum+= exp[j];
+    }
 …
   Free( (ADDRESS) exp, (n+1) * sizeof(Exponent_t) );
   pOrdPolyMerge(pit);
   return pit;
 …
   w= (number *)Alloc( cn * sizeof(number) );
   c= (number *)Alloc( cn * sizeof(number) );
+  for ( j= 0; j < cn; j++ ) {
+  for ( j= 0; j < cn; j++ )
+  {
     w[j]= nInit(0);
     c[j]= nInit(0);
+  }
+  if ( cn == 1 ) {
+  if ( cn == 1 )
+  {
     nDelete( &w[0] );
     w[0]= nCopy(q[0]);
+  } else {
+  }
+  else
+  {
     nDelete( &c[cn-1] );
     c[cn-1]= nCopy(x[0]);
 …
       for ( j= (cn-i-1); j <= (cn-2); j++) { // j=(cn+1-i); j <= (cn-1)
         nDelete( &tmp1 );
         tmp1= nMult( xx, c[j+1] );           // c[j]= c[j] + (xx * c[j+1])
         newnum= nAdd( c[j], tmp1 );
         nDelete( &c[j] );
         c[j]= newnum;
+        nDelete( &tmp1 );
+        tmp1= nMult( xx, c[j+1] );           // c[j]= c[j] + (xx * c[j+1])
+        newnum= nAdd( c[j], tmp1 );
+        nDelete( &c[j] );
+        c[j]= newnum;
+      }
 …
       xx= nCopy(x[i]);                     // xx= x[i]
       nDelete( &t );
+      nDelete( &t );
       t= nInit( 1 );                         // t= b= 1
       nDelete( &b );
       b= nInit( 1 );
+      b= nInit( 1 );
       nDelete( &s );                         // s= q[cn-1]
       s= nCopy( q[cn-1] );
       for ( k= cn-1; k >= 1; k-- ) {         // k=cn; k >= 2
         nDelete( &tmp1 );
         tmp1= nMult( xx, b );                // b= c[k] + (xx * b)
         nDelete( &b );
+        b= nAdd( c[k], tmp1 );
         nDelete( &tmp1 );
         tmp1= nMult( q[k-1], b );            // s= s + (q[k-1] * b)
         newnum= nAdd( s, tmp1 );
         nDelete( &s );
         s= newnum;
         nDelete( &tmp1 );
         tmp1= nMult( xx, t );                // t= (t * xx) + b
         newnum= nAdd( tmp1, b );
         nDelete( &t );
         t= newnum;
+        nDelete( &tmp1 );
+        tmp1= nMult( xx, b );                // b= c[k] + (xx * b)
+        nDelete( &b );
+        b= nAdd( c[k], tmp1 );
+        nDelete( &tmp1 );
+        tmp1= nMult( q[k-1], b );            // s= s + (q[k-1] * b)
+        newnum= nAdd( s, tmp1 );
+        nDelete( &s );
+        s= newnum;
+        nDelete( &tmp1 );
+        tmp1= nMult( xx, t );                // t= (t * xx) + b
+        newnum= nAdd( tmp1, b );
+        nDelete( &t );
+        t= newnum;
+      }
 …
 //-> definitions
 #define MR       8        // never change this value
 #define MT      20
+#define MT      20
 #define MAXIT   (MT*MR)   // max number of iterations in laguer root finder
 …
 //<-
 //-> rootContainer::rootContainer()
+//-> rootContainer::rootContainer()
 rootContainer::rootContainer()
+{
 …
 //<-
 //-> rootContainer::~rootContainer()
+//-> rootContainer::~rootContainer()
 rootContainer::~rootContainer()
+{
 …
   // free coeffs, ievpoint
+  if ( ievpoint != NULL ) {
+  if ( ievpoint != NULL )
+  {
     for ( i=0; i < anz+2; i++ ) nDelete( ievpoint + i );
     Free( (ADDRESS)ievpoint, (anz+2) * sizeof( number ) );
 …
   for ( i=0; i < tdg; i++ ) delete theroots[i];
   Free( (ADDRESS) theroots, (tdg)*sizeof(gmp_complex*) );
   mprPROTnl("~rootContainer()");
+}
 //<-
 //-> void rootContainer::fillContainer( ... )
 void rootContainer::fillContainer( number *_coeffs, number *_ievpoint,
+                                   const int _var, const int _tdg,
                                    const rootType  _rt, const int _anz )
+//-> void rootContainer::fillContainer( ... )
+void rootContainer::fillContainer( number *_coeffs, number *_ievpoint,
+                                   const int _var, const int _tdg,
+                                   const rootType  _rt, const int _anz )
+{
   int i;
 …
   anz=_anz;
+  for ( i=0; i <= tdg; i++ ) {
+    if ( nEqual(coeffs[i],nn) ) {
+  for ( i=0; i <= tdg; i++ )
+  {
+    if ( nEqual(coeffs[i],nn) )
+    {
       nDelete( &coeffs[i] );
       coeffs[i]=NULL;
 …
     ievpoint= (number *)Alloc( (anz+2) * sizeof( number ) );
     for (i=0; i < anz+2; i++) ievpoint[i]= nCopy( _ievpoint[i] );
+  }
+  }
   theroots= NULL;
 …
 //<-
 //-> poly rootContainer::getPoly()
+//-> poly rootContainer::getPoly()
 poly rootContainer::getPoly()
+{
 …
   poly ppos;
+  if ( (rt == cspecial) || ( rt == cspecialmu ) ) {
+    for ( i= tdg; i >= 0; i-- ) {
+      if ( coeffs[i] ) {
+        poly p= pOne();
+        //pSetExp( p, var+1, i);
+        pSetExp( p, 1, i);
+        pSetCoeff( p, nCopy( coeffs[i] ) );
+        pSetm( p );
+        if (result) {
+          ppos->next=p;
+          ppos=ppos->next;
+        } else {
+          result=p;
+          ppos=p;
+        }
+  if ( (rt == cspecial) || ( rt == cspecialmu ) )
+  {
+    for ( i= tdg; i >= 0; i-- )
+    {
+      if ( coeffs[i] )
+      {
+        poly p= pOne();
+        //pSetExp( p, var+1, i);
+        pSetExp( p, 1, i);
+        pSetCoeff( p, nCopy( coeffs[i] ) );
+        pSetm( p );
+        if (result)
+        {
+          ppos->next=p;
+          ppos=ppos->next;
+        }
+        else
+        {
+          result=p;
+          ppos=p;
+        }
+      }
+    }
     pSetm( result );
+  }
+  }
   return result;
 …
 //<-
 //-> const gmp_complex & rootContainer::opterator[] ( const int i )
+//-> const gmp_complex & rootContainer::opterator[] ( const int i )
 // this is now inline!
 //  gmp_complex & rootContainer::operator[] ( const int i )
 //  {
+//    if ( found_roots && ( i >= 0) && ( i < tdg ) ) {
+//    if ( found_roots && ( i >= 0) && ( i < tdg ) )
+//    {
 //      return *theroots[i];
 //    }
+//    }
 //    // warning
 //    Werror("rootContainer::getRoot: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
 …
 //<-
 //-> gmp_complex & rootContainer::evPointCoord( int i )
+//-> gmp_complex & rootContainer::evPointCoord( int i )
 gmp_complex & rootContainer::evPointCoord( const int i )
+{
   if (! ((i >= 0) && (i < anz+2) ) )
     Werror("rootContainer::evPointCoord: index out of range");
+    WerrorS("rootContainer::evPointCoord: index out of range");
   if (ievpoint == NULL)
     Werror("rootContainer::evPointCoord: ievpoint == NULL");
+    WerrorS("rootContainer::evPointCoord: ievpoint == NULL");
   if ( (rt == cspecialmu) && found_roots ) {  // FIX ME
+    if ( ievpoint[i] != NULL ) {
+    if ( ievpoint[i] != NULL )
+    {
       gmp_complex *tmp= new gmp_complex();
       *tmp= numberToGmp_Complex(ievpoint[i]);
       return *tmp;
+    } else {
+    }
+    else
+    {
       Werror("rootContainer::evPointCoord: NULL index %d",i);
+    }
+  }
   // warning
   Werror("rootContainer::evPointCoord: Wrong index %d, found_roots %s",i,found_roots?"true":"false");
 …
 //<-
 //-> bool rootContainer::changeRoots( int from, int to )
+//-> bool rootContainer::changeRoots( int from, int to )
 bool rootContainer::swapRoots( const int from, const int to )
+{
+  if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) ) {
+    if ( to != from ) {
+  if ( found_roots && ( from >= 0) && ( from < tdg ) && ( to >= 0) && ( to < tdg ) )
+  {
+    if ( to != from )
+    {
       gmp_complex tmp( *theroots[from] );
       *theroots[from]= *theroots[to];
 …
+    }
     return true;
+  }
+  }
   // warning
   Werror(" rootContainer::changeRoots: Wrong index %d, %d",from,to);
 …
 //<-
 //-> void rootContainer::solver()
+//-> void rootContainer::solver()
 bool rootContainer::solver( const int polishmode )
+{
 …
   // copy the coefficients of type number to type gmp_complex
   gmp_complex **ad= (gmp_complex**)Alloc( (tdg+1)*sizeof(gmp_complex*) );
+  for ( i=0; i <= tdg; i++ ) {
+  for ( i=0; i <= tdg; i++ )
+  {
     ad[i]= new gmp_complex();
     if ( coeffs[i] ) *ad[i] = numberToGmp_Complex( coeffs[i] );
 …
   // now solve
+  switch (polishmode) {
+  switch (polishmode)
+  {
   case PM_NONE:
   case PM_POLISH:
     found_roots= laguer_driver( ad, theroots, polishmode == PM_POLISH );
+    if (!found_roots) {
+      Werror("rootContainer::solver: No roots found!");
+    if (!found_roots)
+    {
+      WerrorS("rootContainer::solver: No roots found!");
       goto solverend;
+    }
     break;
   case PM_CORRUPT:
     found_roots= laguer_driver( ad, theroots, false );
+    found_roots= laguer_driver( ad, theroots, false );
     // corrupt the roots
     for ( i= 0; i < tdg; i++ )
       *theroots[i]= *theroots[i] * (gmp_float)(1.0+0.01*(mprfloat)i);
     // and interpolate again
+    // and interpolate again
     found_roots= laguer_driver( ad, theroots, true );
+    if (!found_roots) {
+      Werror("rootContainer::solver: No roots found!");
+    if (!found_roots)
+    {
+      WerrorS("rootContainer::solver: No roots found!");
       goto solverend;
+    }
 …
 //<-
 //-> gmp_complex* rootContainer::laguer_driver( bool polish )
+//-> gmp_complex* rootContainer::laguer_driver( bool polish )
 bool rootContainer::laguer_driver(gmp_complex ** a, gmp_complex ** roots, bool polish )
+{
 …
   for ( i=0; i <= tdg; i++ ) ad[i]= new gmp_complex( *a[i] );
+  for ( j= tdg; j >= 1; j-- ) {
+  for ( j= tdg; j >= 1; j-- )
+  {
     x= gmp_complex();
 …
     mprSTICKYPROT(ST_ROOTS_LGSTEP);
     if ( its > MAXIT ) {  // error
       Werror("rootContainer::laguer_driver: To many iterations!");
+      WerrorS("rootContainer::laguer_driver: To many iterations!");
       ret= false;
       goto theend;
+    }
+    if ( abs(x.imag()) <= (gmp_float)(2.0*EPS)*abs(x.real())) {
+    if ( abs(x.imag()) <= (gmp_float)(2.0*EPS)*abs(x.real()))
+    {
       x= gmp_complex( x.real() );
+    }
     *roots[j-1]= x;
     b= *ad[j];
+    for ( jj= j-1; jj >= 0; jj-- ) {
+    for ( jj= j-1; jj >= 0; jj-- )
+    {
       c= *ad[jj];
       *ad[jj]= b;
 …
+  }
+  if ( polish ) {
+  if ( polish )
+  {
     mprSTICKYPROT(ST_ROOTS_LGPOLISH);
     for ( i=0; i <= tdg; i++ ) *ad[i]=*a[i];
+    for ( j= 1; j <= tdg; j++ ) {
+    for ( j= 1; j <= tdg; j++ )
+    {
       // run laguer alg with corrupted roots
       laguer( ad, tdg, roots[j-1], &its );
       mprSTICKYPROT(ST_ROOTS_LGSTEP);
+      if ( its > MAXIT ) {  // error
+        Werror("rootContainer::laguer_driver: To many iterations!");
+        ret= false;
+        goto theend;
+      }
+    }
+    for ( j= 2; j <= tdg; j++ ) {
+      if ( its > MAXIT )
+      {  // error
+        WerrorS("rootContainer::laguer_driver: To many iterations!");
+        ret= false;
+        goto theend;
+      }
+    }
+    for ( j= 2; j <= tdg; j++ )
+    {
       // sort root by their absolute real parts by straight insertion
       x= *roots[j-1];
+      for ( i= j-1; i >= 1; i-- ) {
+        if ( abs(roots[i-1]->real()) <= abs(x.real()) ) break;
+        *roots[i]= *roots[i-1];
+      for ( i= j-1; i >= 1; i-- )
+      {
+        if ( abs(roots[i-1]->real()) <= abs(x.real()) ) break;
+        *roots[i]= *roots[i-1];
+      }
       *roots[i]= x;
 …
 //<-
 //-> void rootContainer::laguer(...)
+//-> void rootContainer::laguer(...)
 void rootContainer::laguer(gmp_complex ** a, int m, gmp_complex *x, int *its)
+{
 …
   gmp_float epss(1.0/pow(10.0,(int)(gmp_output_digits+gmp_output_digits/4)));
+  for ( iter= 1; iter <= MAXIT; iter++ ) {
+  for ( iter= 1; iter <= MAXIT; iter++ )
+  {
     mprSTICKYPROT(ST_ROOTS_LG);
 …
     b= *a[m];
     err_g= abs(b);
+    err_g= abs(b);
     d= gmp_complex();
     f= gmp_complex();
+    abx_g= abs(*x);
+    for ( j= m-1; j >= 0; j-- ) {
+    abx_g= abs(*x);
+    for ( j= m-1; j >= 0; j-- )
+    {
       f= ( *x * f ) + d;
       d= ( *x * d ) + b;
       b= ( *x * b ) + *a[j];
       err_g= abs( b ) + ( abx_g * err_g );
+      err_g= abs( b ) + ( abx_g * err_g );
+    }
     err_g *= epss; // EPSS;
     if ( abs(b) <= err_g ) return;
+    if ( abs(b) <= err_g ) return;
     g= d / b;
 …
     gp= g + sq;
     gm= g - sq;
     abp_g= abs( gp );
     abm_g= abs( gm );
+    abp_g= abs( gp );
+    abm_g= abs( gm );
     if ( abp_g < abm_g ) gp= gm;
     dx = ( (max(abp_g,abm_g) > (gmp_float)0.0)
+           ? ( gmp_complex( (mprfloat)m ) / gp )
+           : ( gmp_complex( cos((mprfloat)iter),sin((mprfloat)iter))
                * exp(log((gmp_float)1.0+abx_g))) );
+    dx = ( (max(abp_g,abm_g) > (gmp_float)0.0)
+           ? ( gmp_complex( (mprfloat)m ) / gp )
+           : ( gmp_complex( cos((mprfloat)iter),sin((mprfloat)iter))
+               * exp(log((gmp_float)1.0+abx_g))) );
     x1= *x - dx;
 …
 //-----------------------------------------------------------------------------
 //-> rootArranger::rootArranger(...)
 rootArranger::rootArranger( rootContainer ** _roots,
+                            rootContainer ** _mu,
                             const int _howclean )
+//-> rootArranger::rootArranger(...)
+rootArranger::rootArranger( rootContainer ** _roots,
+                            rootContainer ** _mu,
+                            const int _howclean )
   : roots(_roots), mu(_mu), howclean(_howclean)
+{
 …
   // find roots of polys given by coeffs in roots
   rc= roots[0]->getAnzElems();
+  for ( i= 0; i < rc; i++ )
+    if ( !roots[i]->solver( howclean ) ) {
+  for ( i= 0; i < rc; i++ )
+    if ( !roots[i]->solver( howclean ) )
+    {
       found_roots= false;
       return;
 …
   // find roots of polys given by coeffs in mu
   mc= mu[0]->getAnzElems();
+  for ( i= 0; i < mc; i++ )
+    if ( ! mu[i]->solver( howclean ) ) {
+  for ( i= 0; i < mc; i++ )
+    if ( ! mu[i]->solver( howclean ) )
+    {
       found_roots= false;
       return;
 …
 //<-
 //-> void rootArranger::arrange()
+//-> void rootArranger::arrange()
 void rootArranger::arrange()
+{
 …
   for ( xkoord= 0; xkoord < anzm; xkoord++ ) {    // für x1,x2, x1,x2,x3, x1,x2,...,xn
     for ( r= 0; r < anzr; r++ ) {                 // für jede Nullstelle
       // (x1-koordinate) * evp[1] + (x2-koordinate) * evp[2] +
+    for ( r= 0; r < anzr; r++ ) {                 // für jede Nullstelle
+      // (x1-koordinate) * evp[1] + (x2-koordinate) * evp[2] +
       //                                  ... + (xkoord-koordinate) * evp[xkoord]
       tmp= gmp_complex();
+      for ( xk =0; xk <= xkoord; xk++ ){
+        tmp -= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); //xk+1
+      for ( xk =0; xk <= xkoord; xk++ )
+      {
+        tmp -= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1); //xk+1
+      }
       found= false;
+      for ( rtest= r; rtest < anzr; rtest++ ) {   // für jede Nullstelle
+        zwerg = tmp - (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xk+1, xkoord+2
+        for ( mtest= 0; mtest < anzr; mtest++ ) {
+          //      if ( tmp == (*mu[xkoord])[mtest] ) {
+          if ( ((zwerg.real() <= (*mu[xkoord])[mtest].real() + mprec) &&
+                (zwerg.real() >= (*mu[xkoord])[mtest].real() - mprec)) &&
+               ((zwerg.imag() <= (*mu[xkoord])[mtest].imag() + mprec) &&
+                (zwerg.imag() >= (*mu[xkoord])[mtest].imag() - mprec)) ) {
+            roots[xk]->swapRoots( r, rtest );
+            found= true;
+            break;
+          }
+        }
+      for ( rtest= r; rtest < anzr; rtest++ ) {   // für jede Nullstelle
+        zwerg = tmp - (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xk+1, xkoord+2
+        for ( mtest= 0; mtest < anzr; mtest++ )
+        {
+          //          if ( tmp == (*mu[xkoord])[mtest] )
+          //          {
+          if ( ((zwerg.real() <= (*mu[xkoord])[mtest].real() + mprec) &&
+                (zwerg.real() >= (*mu[xkoord])[mtest].real() - mprec)) &&
+               ((zwerg.imag() <= (*mu[xkoord])[mtest].imag() + mprec) &&
+                (zwerg.imag() >= (*mu[xkoord])[mtest].imag() - mprec)) )
+           {
+             roots[xk]->swapRoots( r, rtest );
+             found= true;
+             break;
+           }
+        }
       } // rtest
+      if ( !found ) {
+        Werror("rootArranger::arrange: No match? coord %d, root %d.",xkoord,r);
+      if ( !found )
+      {
+        Werror("rootArranger::arrange: No match? coord %d, root %d.",xkoord,r);
 #ifdef mprDEBUG_PROT
+        Werror("One of these ...");
+        for ( rtest= r; rtest < anzr; rtest++ ) {
+          tmp= gmp_complex();
+          for ( xk =0; xk <= xkoord; xk++ ){
+            tmp-= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1);
+          }
+          tmp-= (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xkoord+2
+          Werror("  %s",complexToStr(tmp,gmp_output_digits+1),rtest);
+        }
+        Werror(" ... must match to one of these:");
+        for ( mtest= 0; mtest < anzr; mtest++ ) {
+          Werror("                  %s",complexToStr((*mu[xkoord])[mtest],gmp_output_digits+1));
+        }
+        WerrorS("One of these ...");
+        for ( rtest= r; rtest < anzr; rtest++ )
+        {
+          tmp= gmp_complex();
+          for ( xk =0; xk <= xkoord; xk++ )
+          {
+            tmp-= (*roots[xk])[r] * mu[xkoord]->evPointCoord(xk+1);
+          }
+          tmp-= (*roots[xk])[rtest] * mu[xkoord]->evPointCoord(xk+1); // xkoord+2
+          Werror("  %s",complexToStr(tmp,gmp_output_digits+1),rtest);
+        }
+        WerrorS(" ... must match to one of these:");
+        for ( mtest= 0; mtest < anzr; mtest++ )
+        {
+          Werror("                  %s",complexToStr((*mu[xkoord])[mtest],gmp_output_digits+1));
+        }
 #endif
+      }
 …
 //<-
 //-> lists rootArranger::listOfRoots( int oprec )
+//-> lists rootArranger::listOfRoots( int oprec )
 lists rootArranger::listOfRoots( const unsigned int oprec )
+{
 …
   lists listofroots= (lists)Alloc( sizeof(slists) ); // must be done this way!
+  if ( found_roots ) {
+  if ( found_roots )
+  {
     listofroots->Init( count );
+    for (i=0; i < count; i++) {
+    for (i=0; i < count; i++)
+    {
       lists onepoint= (lists)Alloc(sizeof(slists)); // must be done this way!
       onepoint->Init(elem);
+      for ( j= 0; j < elem; j++ ) {
+        if ( !rField_is_long_C() ) {
+          onepoint->m[j].rtyp=STRING_CMD;
+          onepoint->m[j].data=(void *)complexToStr((*roots[j])[i],oprec);
+        } else {
+          onepoint->m[j].rtyp=NUMBER_CMD;
+          onepoint->m[j].data=(void *)nCopy((number)(roots[j]->getRoot(i)));
+        }
+        onepoint->m[j].next= NULL;
+        onepoint->m[j].name= NULL;
+      for ( j= 0; j < elem; j++ )
+      {
+        if ( !rField_is_long_C() )
+        {
+          onepoint->m[j].rtyp=STRING_CMD;
+          onepoint->m[j].data=(void *)complexToStr((*roots[j])[i],oprec);
+        }
+        else
+        {
+          onepoint->m[j].rtyp=NUMBER_CMD;
+          onepoint->m[j].data=(void *)nCopy((number)(roots[j]->getRoot(i)));
+        }
+        onepoint->m[j].next= NULL;
+        onepoint->m[j].name= NULL;
+      }
       listofroots->m[i].rtyp=LIST_CMD;
 …
+    }
+  } else {
+  }
+  else
+  {
     listofroots->Init( 0 );
+  }
 …
   int i,ip,ir,is,k,kh,kp,m12,nl1,nl2;
   int *l1,*l2,*l3;
+  mprfloat q1, bmax;
+  if ( m != (m1+m2+m3)) {
+  mprfloat q1, bmax;
+  if ( m != (m1+m2+m3))
+  {
     // error: bad input
     error(Werror(" bad input constraint counts in simplex ");)
+    error(WerrorS(" bad input constraint counts in simplex ");)
     *icase=-2;
     return;
 …
   for ( k=1; k<=n; k++ ) l1[k]=izrov[k]=k;
   nl2=m;
+  for ( i=1; i<=m; i++ ) {
+    if ( a[i+1][1] < 0.0 ) {
+  for ( i=1; i<=m; i++ )
+  {
+    if ( a[i+1][1] < 0.0 )
+    {
       // error: bad input
       error(Werror(" bad input tableau in simplex ");)
+      error(WerrorS(" bad input tableau in simplex ");)
       *icase=-2;
       // free mem l1,l2,l3;
 …
   for ( i=1; i<=m2; i++) l3[i]= 1;
   ir= 0;
+  if (m2+m3) {
+  if (m2+m3)
+  {
     ir=1;
+    for ( k=1; k <= (n+1); k++ ) {
+    for ( k=1; k <= (n+1); k++ )
+    {
       q1=0.0;
       for ( i=m1+1; i <= m; i++ ) q1+= a[i+1][k];
 …
+    }
+    do {
+    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
+        // free mem l1,l2,l3;
+        Free( (ADDRESS) l3, (m+1) * sizeof(int) );
+        Free( (ADDRESS) l2, (m+1) * sizeof(int) );
+        Free( (ADDRESS) l1, (n+1) * sizeof(int) );
+        return;
+      } else if ( bmax <= SIMPLEX_EPS && a[m+2][1] <= SIMPLEX_EPS ) {
+        m12= m1+m2+1;
+        if ( m12 <= m ) {
+          for ( ip= m12; ip <= m; ip++ ) {
+            if ( iposv[ip] == (ip+n) ) {
+              simp1(a,ip,l1,nl1,1,&kp,&bmax);
+              if ( fabs(bmax) >= SIMPLEX_EPS)
+                goto one;
+            }
+          }
+        }
+        ir= 0;
+        --m12;
+        if ( m1+1 <= m12 )
+          for ( i=m1+1; i <= m12; i++ )
+            if ( l3[i-m1] == 1 )
+              for ( k=1; k <= n+1; k++ )
+                a[i+1][k] = -a[i+1][k];
+        break;
+      if ( bmax <= SIMPLEX_EPS && a[m+2][1] < -SIMPLEX_EPS )
+      {
+        *icase= -1; // no solution found
+        // free mem l1,l2,l3;
+        Free( (ADDRESS) l3, (m+1) * sizeof(int) );
+        Free( (ADDRESS) l2, (m+1) * sizeof(int) );
+        Free( (ADDRESS) l1, (n+1) * sizeof(int) );
+        return;
+      }
+      else if ( bmax <= SIMPLEX_EPS && a[m+2][1] <= SIMPLEX_EPS )
+      {
+        m12= m1+m2+1;
+        if ( m12 <= m )
+        {
+          for ( ip= m12; ip <= m; ip++ )
+          {
+            if ( iposv[ip] == (ip+n) )
+            {
+              simp1(a,ip,l1,nl1,1,&kp,&bmax);
+              if ( fabs(bmax) >= SIMPLEX_EPS)
+                goto one;
+            }
+          }
+        }
+        ir= 0;
+        --m12;
+        if ( m1+1 <= m12 )
+          for ( i=m1+1; i <= m12; i++ )
+            if ( l3[i-m1] == 1 )
+              for ( k=1; k <= n+1; k++ )
+                a[i+1][k] = -a[i+1][k];
+        break;
+      }
       //#if DEBUG
 …
       //#endif
       simp2(a,n,l2,nl2,&ip,kp,&q1);
+      if ( ip == 0 ) {
+        *icase = -1; // no solution found
+        // free mem l1,l2,l3;
+        Free( (ADDRESS) l3, (m+1) * sizeof(int) );
+        Free( (ADDRESS) l2, (m+1) * sizeof(int) );
+        Free( (ADDRESS) l1, (n+1) * sizeof(int) );
+        return;
+      if ( ip == 0 )
+      {
+        *icase = -1; // no solution found
+        // free mem l1,l2,l3;
+        Free( (ADDRESS) l3, (m+1) * sizeof(int) );
+        Free( (ADDRESS) l2, (m+1) * sizeof(int) );
+        Free( (ADDRESS) l1, (n+1) * sizeof(int) );
+        return;
+      }
     one: simp3(a,m+1,n,ip,kp);
 …
       // print_bmat(a,m+2,n+3);
       // #endif
+      if ( iposv[ip] >= (n+m1+m2+1)) {
+        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];
+      } else {
+        if ( iposv[ip] >= (n+m1+1) ) {
+          kh= iposv[ip]-m1-n;
+          if ( l3[kh] ) {
+            l3[kh]= 0;
+            ++a[m+2][kp+1];
+            for ( i=1; i<= m+2; i++ )
+              a[i][kp+1] = -a[i][kp+1];
+          }
+        }
+      if ( iposv[ip] >= (n+m1+m2+1))
+      {
+        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];
+      }
+      else
+      {
+        if ( iposv[ip] >= (n+m1+1) )
+        {
+          kh= iposv[ip]-m1-n;
+          if ( l3[kh] )
+          {
+            l3[kh]= 0;
+            ++a[m+2][kp+1];
+            for ( i=1; i<= m+2; i++ )
+              a[i][kp+1] = -a[i][kp+1];
+          }
+        }
+      }
       is= izrov[kp];
 …
+  }
   /* end of phase 1, have feasible sol, now optimize it */
+  for (;;) {
+  loop
+  {
     // #if DEBUG
     // print_bmat( a, m+1, n+5);
     // #endif
     simp1(a,0,l1,nl1,0,&kp,&bmax);
+    if (bmax <= /*SIMPLEX_EPS*/0.0) {
+    if (bmax <= /*SIMPLEX_EPS*/0.0)
+    {
       *icase=0; // finite solution found
       // free mem l1,l2,l3
 …
+    }
     simp2(a,n,l2,nl2,&ip,kp,&q1);
+    if (ip == 0) {
+    if (ip == 0)
+    {
       //printf("Unbounded:");
       // #if DEBUG
       //       print_bmat( a, m+1, n+1);
       // #endif
       *icase=1;         /* unbounded */
+      *icase=1;                /* unbounded */
       // free mem
       Free( (ADDRESS) l3, (m+1) * sizeof(int) );
 …
     is= izrov[kp];
     izrov[kp]= iposv[ip];
     iposv[ip]= is;
+    iposv[ip]= is;
   }/*for ;;*/
+}
 …
   mprfloat test;
+  if( nll <= 0) {                       /* init'tion: fixed */
+  if( nll <= 0)
+  {                        /* init'tion: fixed */
     *bmax = 0.0;
     return;
 …
   *kp=ll[1];
   *bmax=a[mm+1][*kp+1];
+  for (k=2;k<=nll;k++) {
+    if (iabf == 0) {
+  for (k=2;k<=nll;k++)
+  {
+    if (iabf == 0)
+    {
       test=a[mm+1][ll[k]+1]-(*bmax);
+      if (test > 0.0) {
+        *bmax=a[mm+1][ll[k]+1];
+        *kp=ll[k];
+      }
+    } else {                    /* abs values: have fixed it */
+      if (test > 0.0)
+      {
+        *bmax=a[mm+1][ll[k]+1];
+        *kp=ll[k];
+      }
+    }
+    else
+    {                        /* abs values: have fixed it */
       test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
+      if (test > 0.0) {
+        *bmax=a[mm+1][ll[k]+1];
+        *kp=ll[k];
+      if (test > 0.0)
+      {
+        *bmax=a[mm+1][ll[k]+1];
+        *kp=ll[k];
+      }
+    }
 …
   *ip= 0;
+  for ( i=1; i <= nl2; i++ ) {
+    if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS ) {
+  for ( i=1; i <= nl2; i++ )
+  {
+    if ( a[l2[i]+1][kp+1] < -SIMPLEX_EPS )
+    {
       *q1= -a[l2[i]+1][1] / a[l2[i]+1][kp+1];
       *ip= l2[i];
+      for ( i= i+1; i <= nl2; i++ ) {
+        ii= l2[i];
+        if (a[ii+1][kp+1] < -SIMPLEX_EPS) {
+          q= -a[ii+1][1] / a[ii+1][kp+1];
+          if (q - *q1 < -SIMPLEX_EPS) {
+            *ip=ii;
+            *q1=q;
+          } else if (q - *q1 < SIMPLEX_EPS) {
+            for ( k=1; k<= n; k++ ) {
+              qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
+              q0= -a[ii+1][k+1]/a[ii+1][kp+1];
+              if ( q0 != qp ) break;
+            }
+            if ( q0 < qp ) *ip= ii;
+          }
+        }
+      for ( i= i+1; i <= nl2; i++ )
+      {
+        ii= l2[i];
+        if (a[ii+1][kp+1] < -SIMPLEX_EPS)
+        {
+          q= -a[ii+1][1] / a[ii+1][kp+1];
+          if (q - *q1 < -SIMPLEX_EPS)
+          {
+            *ip=ii;
+            *q1=q;
+          }
+          else if (q - *q1 < SIMPLEX_EPS)
+          {
+            for ( k=1; k<= n; k++ )
+            {
+              qp= -a[*ip+1][k+1]/a[*ip+1][kp+1];
+              q0= -a[ii+1][k+1]/a[ii+1][kp+1];
+              if ( q0 != qp ) break;
+            }
+            if ( q0 < qp ) *ip= ii;
+          }
+        }
+      }
+    }
 …
   piv= 1.0 / a[ip+1][kp+1];
+  for ( ii=1; ii <= i1+1; ii++ )
+    if ( ii -1 != ip ) {
+  for ( ii=1; ii <= i1+1; ii++ )
+  {
+    if ( ii -1 != ip )
+    {
       a[ii][kp+1] *= piv;
+      for ( kk=1; kk <= k1+1; kk++ )
+        if ( kk-1 != kp )
+          a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
+    }
+  for ( kk=1; kk<= k1+1; kk++ )
+      for ( kk=1; kk <= k1+1; kk++ )
+        if ( kk-1 != kp )
+          a[ii][kk] -= a[ip+1][kk] * a[ii][kp+1];
+    }
+  }
+  for ( kk=1; kk<= k1+1; kk++ )
     if ( kk-1 != kp ) a[ip+1][kk] *= -piv;
   a[ip+1][kp+1]= piv;
 …
 // compile-command-2: "make install" ***
 // End: ***

Singular/mpr_numeric.h

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: mpr_numeric.h,v 1.1 1999-06-28 12:48:16 wenk Exp $ */
+/* $Id: mpr_numeric.h,v 1.2 1999-06-28 16:06:30 Singular Exp $ */
 /*
+/*
 * ABSTRACT - multipolynomial resultants - numeric stuff
 *            ( root finder, vandermonde system solver, simplex )
+*            ( root finder, vandermonde system solver, simplex )
 */
 …
 // define polish mode when finding roots
 #define PM_NONE    1
+#define PM_NONE    1
 #define PM_POLISH  2
 #define PM_CORRUPT 3
 …
+{
 public:
   vandermonde( const long _cn, const long _n,
+  vandermonde( const long _cn, const long _n,
                const long _maxdeg, number *_p, const bool _homog = true );
   ~vandermonde();
   /** Solves the Vandermode linear system
    *    \sum_{i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
+  /** Solves the Vandermode linear system
+   *    \sum_{i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
      * Any computations are done using type number to get high pecision results.
    * @param  q n-tuple of results (right hand of equations)
 …
   ~rootContainer();
   void fillContainer( number *_coeffs, number *_ievpoint,
                       const int _var, const int _tdg,
+  void fillContainer( number *_coeffs, number *_ievpoint,
+                      const int _var, const int _tdg,
                       const rootType _rt, const int _anz );
 …
   rootContainer( const rootContainer & v );
   /** Given the degree tdg and the tdg+1 complex coefficients ad[0..tdg]
+  /** Given the degree tdg and the tdg+1 complex coefficients ad[0..tdg]
    * (generated from the number coefficients coeffs[0..tdg]) of the polynomial
    * this routine successively calls "laguer" and finds all m complex roots in
+   * this routine successively calls "laguer" and finds all m complex roots in
    * roots[0..tdg]. The bool var "polish" should be input as "true" if polishing
    * (also by "laguer") is desired, "false" if the roots will be subsequently
 …
    */
   void laguer(gmp_complex ** a, int m, gmp_complex * x, int * its);
   int var;
   int tdg;
 …
   gmp_complex ** theroots;
   int anz;
   bool found_roots;
 …
 //-> class rootArranger
 class rootArranger
+class rootArranger
+{
 public:
   rootArranger( rootContainer ** _roots,
                 rootContainer ** _mu,
+public:
+  rootArranger( rootContainer ** _roots,
+                rootContainer ** _mu,
                 const int _howclean = PM_CORRUPT );
   ~rootArranger() {}
 …
 // compile-command-1: "make installg" ***
 // compile-command-2: "make install" ***
+// End: ***
+// End: ***

Singular/pcv.cc

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR      *
 *****************************************/
 /* $Id: pcv.cc,v 1.25 1999-06-11 17:13:40 mschulze Exp $ */
+/* $Id: pcv.cc,v 1.26 1999-06-28 16:06:30 Singular Exp $ */
 /*
 * ABSTRACT: conversion between polys and coef vectors
 …
       l0->m[i].data=pCopy((poly)l1->m[i].data);
       if(i<=l2->nr&&l2->m[i].rtyp==l1->m[i].rtyp)
         l0->m[i].data=pAdd(l0->m[i].data,pCopy((poly)l2->m[i].data));
+        l0->m[i].data=pAdd((poly)l0->m[i].data,pCopy((poly)l2->m[i].data));
+    }
     else
 …
+  {
     k=j;
     for(j=0;j<pcvMaxDegree&&pcvIndex[i][j]<=n;j++);
+    for(j=0; (j<pcvMaxDegree) && (pcvIndex[i][j]<=(unsigned)n); j++);
     j--;
     n-=pcvIndex[i][j];

Singular/shortfl.cc

-                      r0f5091
+                      re858e7
 *  Computer Algebra System SINGULAR     *
 ****************************************/
 /* $Id: shortfl.cc,v 1.8 1999-05-10 15:10:55 Singular Exp $ */
+/* $Id: shortfl.cc,v 1.9 1999-06-28 16:06:31 Singular Exp $ */
 /*
 …
     if(i>4)
+    {
       Werror("float overflow");
+      WerrorS("float overflow");
       return nf(rr).N();
+    }
 …
+  {
     if(j==s)
       Werror("float overflow");
+      WerrorS("float overflow");
     return nf(rr).N();
+  }
 …
       MPZ_CLEAR(g);
       if(j==s)
         Werror("float overflow");
+        WerrorS("float overflow");
       return nf(rr).N();
+    }

Context Navigation

Changeset e858e7 in git

Legend:

Singular/clapsing.cc

Singular/intvec.cc

Singular/intvec.h

Singular/mpr_base.cc

Singular/mpr_base.h

Singular/mpr_complex.cc

Singular/mpr_complex.h

Singular/mpr_global.h

Singular/mpr_inout.cc

Singular/mpr_inout.h

Singular/mpr_numeric.cc

Singular/mpr_numeric.h

Singular/pcv.cc

Singular/shortfl.cc

Download in other formats: