Changeset e9ade0f in git

Timestamp:

Feb 22, 2007, 3:46:51 PM (16 years ago)

Author:

Michael Brickenstein <bricken@…>

Branches:

(u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')

Children:

05023393670050d3a3678f07351eda09ade7a6eb

Parents:

ebe987f7339e431a23e717106ffe22dacf14a1a2

Message:

*bricken: templates for matrix


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

Files:

• Singular/claptmpl.cc (modified) (2 diffs)
• kernel/tgb.cc (modified) (9 diffs)
• kernel/tgb_internal.h (modified) (2 diffs)

Singular/claptmpl.cc

@@ -3 +3 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-// $Id: claptmpl.cc,v 1.39 2007-02-22 10:45:12 bricken Exp $
+// $Id: claptmpl.cc,v 1.40 2007-02-22 14:46:51 bricken Exp $
 /*
 * ABSTRACT - instantiation of all templates
@@ -227 +227 @@
 //template class std::vector<std::vector<NoroPlaceHolder> >;
 template class std::vector<DataNoroCacheNode*>;
-template SparseRow* noro_red_to_non_poly_t<unsigned char>(poly p, int &len, NoroCache* cache,slimgb_alg* c);
-template SparseRow* noro_red_to_non_poly_t<unsigned short>(poly p, int &len, NoroCache* cache,slimgb_alg* c);
-template SparseRow* noro_red_to_non_poly_t<unsigned int>(poly p, int &len, NoroCache* cache,slimgb_alg* c);
+template SparseRow* noro_red_to_non_poly_t<tgb_uint16>(poly p, int &len, NoroCache* cache,slimgb_alg* c);
+template SparseRow* noro_red_to_non_poly_t<tgb_uint32>(poly p, int &len, NoroCache* cache,slimgb_alg* c);
+template SparseRow* noro_red_to_non_poly_t<tgb_uint8>(poly p, int &len, NoroCache* cache,slimgb_alg* c);
+template void simplest_gauss_modp<tgb_uint16> (tgb_uint16* a, int nrows,int ncols);
+template void simplest_gauss_modp<tgb_uint32> (tgb_uint32* a, int nrows,int ncols);
+template void simplest_gauss_modp<tgb_uint8> (tgb_uint8* a, int nrows,int ncols);
+template poly row_to_poly<tgb_uint8>(tgb_uint8* row, poly* terms, int tn, ring r);
+template poly row_to_poly<tgb_uint32>(tgb_uint32* row, poly* terms, int tn, ring r);
+template poly row_to_poly<tgb_uint16>(tgb_uint16* row, poly* terms, int tn, ring r);
 //std::priority_queue<MonRedRes>
 #endif

kernel/tgb.cc

@@ -5 +5 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: tgb.cc,v 1.146 2007-02-22 13:56:32 bricken Exp $ */
+/* $Id: tgb.cc,v 1.147 2007-02-22 14:46:49 bricken Exp $ */
 /*
 * ABSTRACT: slimgb and F4 implementation
@@ -1827 +1827 @@
   return TRUE;
 }
-static int terms_sort_crit(const void* a, const void* b){
+int terms_sort_crit(const void* a, const void* b){
   return -pLmCmp(*((poly*) a),*((poly*) b));
 }
@@ -1866 +1866 @@
 }
 typedef unsigned short number_type;
-#define prec_cast(a) (unsigned int) (a)
-#define to_number_type(a) (number_type) prec_cast(a)
-int modP_lastIndexRow(number_type* row,int ncols){
-  int lastIndex;
-  const number_type zero=0;//npInit(0);
-  for(lastIndex=ncols-1;lastIndex>=0;lastIndex--){
-    if (!(row[lastIndex]==zero)){
-      return lastIndex;
-    }
-  }
-  return -1;
-}
-
-class ModPMatrixProxyOnArray{
-public:
-  friend class ModPMatrixBackSubstProxOnArray;
-  int ncols,nrows;
-  ModPMatrixProxyOnArray(number_type* array, int nrows, int ncols){
-    this->ncols=ncols;
-    this->nrows=nrows;
-    rows=(number_type**) omalloc(nrows*sizeof(number_type*));
-    startIndices=(int*)omalloc(nrows*sizeof(int));
-    int i;
-    for(i=0;i<nrows;i++){
-      rows[i]=array+(i*ncols);
-      updateStartIndex(i,-1);
-    }
-  }
-  ~ModPMatrixProxyOnArray(){
-    omfree(rows);
-    omfree(startIndices);
-  }
-  
-  void permRows(int i, int j){
-    number_type* h=rows[i];
-    rows[i]=rows[j];
-    rows[j]=h;
-    int hs=startIndices[i];
-    startIndices[i]=startIndices[j];
-    startIndices[j]=hs;
-  }
-  void multiplyRow(int row, number_type coef){
-    int i;
-    number_type* row_array=rows[row];
-    for(i=startIndices[row];i<ncols;i++){
-      row_array[i]=to_number_type(npMult((number) row_array[i],(number) coef));
-    }
-  }
-  void reduceOtherRowsForward(int r){
-
-    //assume rows "under r" have bigger or equal start index
-    number_type* row_array=rows[r];
-    number_type zero=to_number_type(npInit(0));
-    int start=startIndices[r];
-    number_type coef=row_array[start];
-    assume(start<ncols);
-    int other_row;
-    assume(!(npIsZero((number) row_array[start])));
-    if (!(npIsOne((number) coef)))
-      multiplyRow(r,to_number_type(npInvers((number) coef)));
-    assume(npIsOne((number) row_array[start]));
-    int lastIndex=modP_lastIndexRow(row_array, ncols);
-    number minus_one=npInit(-1);
-    for (other_row=r+1;other_row<nrows;other_row++){
-      assume(startIndices[other_row]>=start);
-      if (startIndices[other_row]==start){
-        int i;
-        number_type* other_row_array=rows[other_row];
-        number coef2=npNeg((number) other_row_array[start]);
-        if (coef2==minus_one){
-          for(i=start;i<=lastIndex;i++){
-            if (row_array[i]!=zero)
-              other_row_array[i]=to_number_type(npSubM((number) other_row_array[i], (number) row_array[i]));
-          }
-      }else {
-          //assume(FALSE);
-          for(i=start;i<=lastIndex;i++){
-            if (row_array[i]!=zero)
-            other_row_array[i]=to_number_type(npAddM(npMult(coef2,(number) row_array[i]),(number) other_row_array[i]));
-          }
-        }
-        updateStartIndex(other_row,start);
-        assume(npIsZero((number) other_row_array[start]));
-      }
-    }
-  }
-  void updateStartIndex(int row,int lower_bound){
-    number_type* row_array=rows[row];
-    assume((lower_bound<0)||(npIsZero((number) row_array[lower_bound])));
-    int i;
-    //number_type zero=npInit(0);
-    for(i=lower_bound+1;i<ncols;i++){
-      if (!(row_array[i]==0))
-        break;
-    }
-    startIndices[row]=i;
-  }
-  int getStartIndex(int row){
-    return startIndices[row];
-  }
-  BOOLEAN findPivot(int &r, int &c){
-    //row>=r, col>=c
-    
-    while(c<ncols){
-      int i;
-      for(i=r;i<nrows;i++){
-        assume(startIndices[i]>=c);
-        if (startIndices[i]==c){
-          //r=i;
-          if (r!=i)
-            permRows(r,i);
-          return TRUE;
-        }
-      }
-      c++;
-    }
-    return FALSE;
-  }
-protected:
-  number_type** rows;
-  int* startIndices;
-};
-class ModPMatrixBackSubstProxOnArray{
-  int *startIndices;
-  number_type** rows;
-  int *lastReducibleIndices;
-  int ncols;
-  int nrows;
-  int nonZeroUntil;
-public:
-  void multiplyRow(int row, number_type coef){
-    int i;
-    number_type* row_array=rows[row];
-    for(i=startIndices[row];i<ncols;i++){
-      row_array[i]=to_number_type(npMult((number) row_array[i],(number) coef));
-    }
-  }
-  ModPMatrixBackSubstProxOnArray(ModPMatrixProxyOnArray& p){
-//  (number_type* array, int nrows, int ncols, int* startIndices, number_type** rows){
-    //we borrow some parameters ;-)
-    //we assume, that nobody changes the order of the rows
-    this->startIndices=p.startIndices;
-    this->rows=p.rows;
-    this->ncols=p.ncols;
-    this->nrows=p.nrows;
-    lastReducibleIndices=(int*) omalloc(nrows*sizeof(int));
-    nonZeroUntil=0;
-    while(nonZeroUntil<nrows){
-      if (startIndices[nonZeroUntil]<ncols){
-       
-        nonZeroUntil++;
-      } else break;
-      
-    }
-    if (TEST_OPT_PROT)
-      Print("rank:%i\n",nonZeroUntil);
-    nonZeroUntil--;
-    int i;
-    for(i=0;i<=nonZeroUntil;i++){
-      assume(startIndices[i]<ncols);
-      assume(!(npIsZero((number) rows[i][startIndices[i]])));
-      assume(startIndices[i]>=i);
-      updateLastReducibleIndex(i,nonZeroUntil+1);
-    }
-  }
-  void updateLastReducibleIndex(int r, int upper_bound){
-    number_type* row_array=rows[r];
-    if (upper_bound>nonZeroUntil) upper_bound=nonZeroUntil+1;
-    int i;
-    const number_type zero=0;//npInit(0);
-    for(i=upper_bound-1;i>r;i--){
-      int start=startIndices[i];
-      assume(start<ncols);
-      if (!(row_array[start]==zero)){
-        lastReducibleIndices[r]=start;
-        return;
-      }
-    }
-    lastReducibleIndices[r]=-1;
-  }
-  void backwardSubstitute(int r){
-    int start=startIndices[r];
-    assume(start<ncols);
-    number_type zero=0;//npInit(0);
-    number_type* row_array=rows[r];
-    assume((!(npIsZero((number) row_array[start]))));
-    assume(start<ncols);
-    int other_row;
-    if (!(npIsOne((number) row_array[r]))){
-      //it should be one, but this safety is not expensive
-      multiplyRow(r, to_number_type(npInvers((number) row_array[start])));
-    }
-    int lastIndex=modP_lastIndexRow(row_array, ncols);
-    assume(lastIndex<ncols);
-    assume(lastIndex>=0);
-    for(other_row=r-1;other_row>=0;other_row--){
-      assume(lastReducibleIndices[other_row]<=start);
-      if (lastReducibleIndices[other_row]==start){
-        number_type* other_row_array=rows[other_row];
-        number coef=npNeg((number) other_row_array[start]);
-        assume(!(npIsZero(coef)));
-        int i;
-        assume(start>startIndices[other_row]);
-        for(i=start;i<=lastIndex;i++){
-          if (row_array[i]!=zero)
-            other_row_array[i]=to_number_type(npAddM(npMult(coef,(number)row_array[i]),(number)other_row_array[i]));
-        }
-        updateLastReducibleIndex(other_row,r);
-      }
-    }
-  }
-  ~ModPMatrixBackSubstProxOnArray(){
-    omfree(lastReducibleIndices);
-  }
-  void backwardSubstitute(){
-    int i;
-    for(i=nonZeroUntil;i>0;i--){
-      backwardSubstitute(i);
-    }
-  }
-};
-void simplest_gauss_modp(number_type* a, int nrows,int ncols){
-  //use memmoves for changing rows
-  if (TEST_OPT_PROT)
-    PrintS("StartGauss\n");
-  ModPMatrixProxyOnArray mat(a,nrows,ncols);
-  
-  int c=0;
-  int r=0;
-  while(mat.findPivot(r,c)){
-    //int pivot=find_pivot()
-      mat.reduceOtherRowsForward(r);
-    r++;
-    c++;
-  }
-  ModPMatrixBackSubstProxOnArray backmat(mat);
-  backmat.backwardSubstitute();
-  //backward substitutions
-  if (TEST_OPT_PROT)
-    PrintS("StopGauss\n");
-}
-static void write_poly_to_row(number_type* row, poly h, poly*terms, int tn, ring r){
-  //poly* base=row;
-  while(h!=NULL){
-    //Print("h:%i\n",h);
-    number coef=p_GetCoeff(h,r);
-    poly* ptr_to_h=(poly*) bsearch(&h,terms,tn,sizeof(poly),terms_sort_crit);
-    assume(ptr_to_h!=NULL);
-    int pos=ptr_to_h-terms;
-    row[pos]=to_number_type(coef);
-    //number_type_array[base+pos]=coef;
-    pIter(h);
-  }
-}
-static poly row_to_poly(number_type* row, poly* terms, int tn, ring r){
-  poly h=NULL;
-  int j;
-  number_type zero=0;//;npInit(0);
-  for(j=tn-1;j>=0;j--){
-    if (!(zero==(row[j]))){
-      poly t=terms[j];
-      t=p_LmInit(t,r);
-      p_SetCoeff(t,(number) row[j],r);
-      pNext(t)=h;
-      h=t;
-    }
-    
-  }
-  return h;
-}
+
+
 #ifdef USE_NORO
 #ifndef NORO_CACHE
@@ -2781 +2513 @@
   return -pLmCmp(((TermNoroDataNode*) a)->t,((TermNoroDataNode*) b)->t);
 }
+
 void noro_step(poly*p,int &pn,slimgb_alg* c){
   //Print("Input rows %d\n",pn);
@@ -2789 +2522 @@
 
   NoroCache cache;
-#ifndef NORO_NON_POLY
-  std::vector<std::vector<NoroPlaceHolder> > place_holders(pn);
-#else
+
   SparseRow** srows=(SparseRow**) omalloc(pn*sizeof(SparseRow*));
-#endif
+
   for(j=0;j<pn;j++){
    
@@ -2801 +2532 @@
     //number coef;
 
-    #ifndef NORO_NON_POLY
-    std::vector<NoroPlaceHolder> ph=noro_red(h,h_len,&cache,c);
-    place_holders[j]=ph;
-    #else
+
     srows[j]=noro_red_to_non_poly(h,h_len,&cache,c);
-    #endif
+
   }
   std::vector<DataNoroCacheNode*> irr_nodes;
@@ -2831 +2559 @@
   qsort(term_nodes,n,sizeof(TermNoroDataNode),term_nodes_sort_crit);
   poly* terms=(poly*) omalloc(n*sizeof(poly));
-  #ifndef NORO_NON_POLY
-  for(j=0;j<n;j++){
-    term_nodes[j].node->term_index=j;
-    terms[j]=term_nodes[j].t;
-  }
-  #else
+
   int* old_to_new_indices=(int*) omalloc(cache.nIrreducibleMonomials*sizeof(int));
   for(j=0;j<n;j++){
@@ -2843 +2566 @@
     terms[j]=term_nodes[j].t;
   }
-  #endif
+
   //if (TEST_OPT_PROT)
   //  Print("Evaluate Rows \n");
-  #ifndef NORO_NON_POLY
-  cache.evaluateRows();
-  #endif
+
   number_type* number_array=(number_type*) omalloc(n*pn*sizeof(number_type));
   memset(number_array,0,sizeof(number_type)*n*pn);
   number zero=npInit(0);
-  if (TEST_OPT_PROT)
-     Print("Evaluate Place Holders\n");
+
   for(j=0;j<pn;j++){
     int i;
@@ -2860 +2580 @@
       row[i]=zero;
     }*/
-    #ifndef NORO_NON_POLY
-    cache.evaluatePlaceHolder(row,place_holders[j]);
-    #else
+
     SparseRow* srow=srows[j];
     for(i=0;i<srow->len;i++){
       int idx=old_to_new_indices[srow->idx_array[i]];
-      row[idx]=to_number_type(srow->coef_array[i]);
+      row[idx]=F4mat_to_number_type(srow->coef_array[i]);
     }
     delete srow;
-    #endif
+
   }
   

kernel/tgb_internal.h

@@ -5 +5 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: tgb_internal.h,v 1.56 2007-02-22 10:45:00 bricken Exp $ */
+/* $Id: tgb_internal.h,v 1.57 2007-02-22 14:46:50 bricken Exp $ */
 /*
  * ABSTRACT: tgb internal .h file
@@ -739 +739 @@
 static wlen_type pair_weighted_length(int i, int j, slimgb_alg* c);
 wlen_type pELength(poly p, ring r);
-void simplest_gauss_modp(number* a, int nrows,int ncols);
-
+int terms_sort_crit(const void* a, const void* b);
+//void simplest_gauss_modp(number* a, int nrows,int ncols);
 // a: a[0,0],a[0,1]....a[nrows-1,ncols-1]
 // assume: field is Zp
-#endif
+#ifdef USE_NORO
+#define slim_prec_cast(a) (unsigned int) (a)
+#define F4mat_to_number_type(a) (number_type) slim_prec_cast(a)
+
+template <class number_type > void write_poly_to_row(number_type* row, poly h, poly*terms, int tn, ring r){
+  //poly* base=row;
+  while(h!=NULL){
+    //Print("h:%i\n",h);
+    number coef=p_GetCoeff(h,r);
+    poly* ptr_to_h=(poly*) bsearch(&h,terms,tn,sizeof(poly),terms_sort_crit);
+    assume(ptr_to_h!=NULL);
+    int pos=ptr_to_h-terms;
+    row[pos]=F4mat_to_number_type(coef);
+    //number_type_array[base+pos]=coef;
+    pIter(h);
+  }
+}
+template <class number_type > poly row_to_poly(number_type* row, poly* terms, int tn, ring r){
+  poly h=NULL;
+  int j;
+  number_type zero=0;//;npInit(0);
+  for(j=tn-1;j>=0;j--){
+    if (!(zero==(row[j]))){
+      poly t=terms[j];
+      t=p_LmInit(t,r);
+      p_SetCoeff(t,(number) row[j],r);
+      pNext(t)=h;
+      h=t;
+    }
+    
+  }
+  return h;
+}
+template <class number_type > int modP_lastIndexRow(number_type* row,int ncols){
+  int lastIndex;
+  const number_type zero=0;//npInit(0);
+  for(lastIndex=ncols-1;lastIndex>=0;lastIndex--){
+    if (!(row[lastIndex]==zero)){
+      return lastIndex;
+    }
+  }
+  return -1;
+}
+template <class number_type>class ModPMatrixBackSubstProxyOnArray;
+template <class number_type > class ModPMatrixProxyOnArray{
+public:
+  friend class ModPMatrixBackSubstProxyOnArray<number_type>;
+               
+  int ncols,nrows;
+  ModPMatrixProxyOnArray(number_type* array, int nrows, int ncols){
+    this->ncols=ncols;
+    this->nrows=nrows;
+    rows=(number_type**) omalloc(nrows*sizeof(number_type*));
+    startIndices=(int*)omalloc(nrows*sizeof(int));
+    int i;
+    for(i=0;i<nrows;i++){
+      rows[i]=array+(i*ncols);
+      updateStartIndex(i,-1);
+    }
+  }
+  ~ModPMatrixProxyOnArray(){
+    omfree(rows);
+    omfree(startIndices);
+  }
+  
+  void permRows(int i, int j){
+    number_type* h=rows[i];
+    rows[i]=rows[j];
+    rows[j]=h;
+    int hs=startIndices[i];
+    startIndices[i]=startIndices[j];
+    startIndices[j]=hs;
+  }
+  void multiplyRow(int row, number_type coef){
+    int i;
+    number_type* row_array=rows[row];
+    for(i=startIndices[row];i<ncols;i++){
+      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
+    }
+  }
+  void reduceOtherRowsForward(int r){
+
+    //assume rows "under r" have bigger or equal start index
+    number_type* row_array=rows[r];
+    number_type zero=F4mat_to_number_type(npInit(0));
+    int start=startIndices[r];
+    number_type coef=row_array[start];
+    assume(start<ncols);
+    int other_row;
+    assume(!(npIsZero((number) row_array[start])));
+    if (!(npIsOne((number) coef)))
+      multiplyRow(r,F4mat_to_number_type(npInvers((number) coef)));
+    assume(npIsOne((number) row_array[start]));
+    int lastIndex=modP_lastIndexRow(row_array, ncols);
+    number minus_one=npInit(-1);
+    for (other_row=r+1;other_row<nrows;other_row++){
+      assume(startIndices[other_row]>=start);
+      if (startIndices[other_row]==start){
+        int i;
+        number_type* other_row_array=rows[other_row];
+        number coef2=npNeg((number) other_row_array[start]);
+        if (coef2==minus_one){
+          for(i=start;i<=lastIndex;i++){
+            if (row_array[i]!=zero)
+              other_row_array[i]=F4mat_to_number_type(npSubM((number) other_row_array[i], (number) row_array[i]));
+          }
+      }else {
+          //assume(FALSE);
+          for(i=start;i<=lastIndex;i++){
+            if (row_array[i]!=zero)
+            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef2,(number) row_array[i]),(number) other_row_array[i]));
+          }
+        }
+        updateStartIndex(other_row,start);
+        assume(npIsZero((number) other_row_array[start]));
+      }
+    }
+  }
+  void updateStartIndex(int row,int lower_bound){
+    number_type* row_array=rows[row];
+    assume((lower_bound<0)||(npIsZero((number) row_array[lower_bound])));
+    int i;
+    //number_type zero=npInit(0);
+    for(i=lower_bound+1;i<ncols;i++){
+      if (!(row_array[i]==0))
+        break;
+    }
+    startIndices[row]=i;
+  }
+  int getStartIndex(int row){
+    return startIndices[row];
+  }
+  BOOLEAN findPivot(int &r, int &c){
+    //row>=r, col>=c
+    
+    while(c<ncols){
+      int i;
+      for(i=r;i<nrows;i++){
+        assume(startIndices[i]>=c);
+        if (startIndices[i]==c){
+          //r=i;
+          if (r!=i)
+            permRows(r,i);
+          return TRUE;
+        }
+      }
+      c++;
+    }
+    return FALSE;
+  }
+protected:
+  number_type** rows;
+  int* startIndices;
+};
+template <class number_type > class ModPMatrixBackSubstProxyOnArray{
+  int *startIndices;
+  number_type** rows;
+  int *lastReducibleIndices;
+  int ncols;
+  int nrows;
+  int nonZeroUntil;
+public:
+  void multiplyRow(int row, number_type coef){
+    int i;
+    number_type* row_array=rows[row];
+    for(i=startIndices[row];i<ncols;i++){
+      row_array[i]=F4mat_to_number_type(npMult((number) row_array[i],(number) coef));
+    }
+  }
+  ModPMatrixBackSubstProxyOnArray<number_type> (ModPMatrixProxyOnArray<number_type> & p){
+//  (number_type* array, int nrows, int ncols, int* startIndices, number_type** rows){
+    //we borrow some parameters ;-)
+    //we assume, that nobody changes the order of the rows
+    this->startIndices=p.startIndices;
+    this->rows=p.rows;
+    this->ncols=p.ncols;
+    this->nrows=p.nrows;
+    lastReducibleIndices=(int*) omalloc(nrows*sizeof(int));
+    nonZeroUntil=0;
+    while(nonZeroUntil<nrows){
+      if (startIndices[nonZeroUntil]<ncols){
+       
+        nonZeroUntil++;
+      } else break;
+      
+    }
+    if (TEST_OPT_PROT)
+      Print("rank:%i\n",nonZeroUntil);
+    nonZeroUntil--;
+    int i;
+    for(i=0;i<=nonZeroUntil;i++){
+      assume(startIndices[i]<ncols);
+      assume(!(npIsZero((number) rows[i][startIndices[i]])));
+      assume(startIndices[i]>=i);
+      updateLastReducibleIndex(i,nonZeroUntil+1);
+    }
+  }
+  void updateLastReducibleIndex(int r, int upper_bound){
+    number_type* row_array=rows[r];
+    if (upper_bound>nonZeroUntil) upper_bound=nonZeroUntil+1;
+    int i;
+    const number_type zero=0;//npInit(0);
+    for(i=upper_bound-1;i>r;i--){
+      int start=startIndices[i];
+      assume(start<ncols);
+      if (!(row_array[start]==zero)){
+        lastReducibleIndices[r]=start;
+        return;
+      }
+    }
+    lastReducibleIndices[r]=-1;
+  }
+  void backwardSubstitute(int r){
+    int start=startIndices[r];
+    assume(start<ncols);
+    number_type zero=0;//npInit(0);
+    number_type* row_array=rows[r];
+    assume((!(npIsZero((number) row_array[start]))));
+    assume(start<ncols);
+    int other_row;
+    if (!(npIsOne((number) row_array[r]))){
+      //it should be one, but this safety is not expensive
+      multiplyRow(r, F4mat_to_number_type(npInvers((number) row_array[start])));
+    }
+    int lastIndex=modP_lastIndexRow(row_array, ncols);
+    assume(lastIndex<ncols);
+    assume(lastIndex>=0);
+    for(other_row=r-1;other_row>=0;other_row--){
+      assume(lastReducibleIndices[other_row]<=start);
+      if (lastReducibleIndices[other_row]==start){
+        number_type* other_row_array=rows[other_row];
+        number coef=npNeg((number) other_row_array[start]);
+        assume(!(npIsZero(coef)));
+        int i;
+        assume(start>startIndices[other_row]);
+        for(i=start;i<=lastIndex;i++){
+          if (row_array[i]!=zero)
+            other_row_array[i]=F4mat_to_number_type(npAddM(npMult(coef,(number)row_array[i]),(number)other_row_array[i]));
+        }
+        updateLastReducibleIndex(other_row,r);
+      }
+    }
+  }
+  ~ModPMatrixBackSubstProxyOnArray<number_type>(){
+    omfree(lastReducibleIndices);
+  }
+  void backwardSubstitute(){
+    int i;
+    for(i=nonZeroUntil;i>0;i--){
+      backwardSubstitute(i);
+    }
+  }
+};
+template <class number_type > void simplest_gauss_modp(number_type* a, int nrows,int ncols){
+  //use memmoves for changing rows
+  if (TEST_OPT_PROT)
+    PrintS("StartGauss\n");
+  ModPMatrixProxyOnArray<number_type> mat(a,nrows,ncols);
+  
+  int c=0;
+  int r=0;
+  while(mat.findPivot(r,c)){
+    //int pivot=find_pivot()
+      mat.reduceOtherRowsForward(r);
+    r++;
+    c++;
+  }
+  ModPMatrixBackSubstProxyOnArray<number_type> backmat(mat);
+  backmat.backwardSubstitute();
+  //backward substitutions
+  if (TEST_OPT_PROT)
+    PrintS("StopGauss\n");
+}
+#endif
+
+#endif