diff --git a/build/pkgs/meataxe/patches/UseErrorPropagation.patch b/build/pkgs/meataxe/patches/UseErrorPropagation.patch
new file mode 100644
index 00000000000..00745e8b967
--- /dev/null
+++ b/build/pkgs/meataxe/patches/UseErrorPropagation.patch
@@ -0,0 +1,1147 @@
+In functions that appear in matrix arithmetic, use specific return values
+on error, and propagate errors. This is *not* done in other parts of
+MeatAxe (e.g., not for greased matrices or polynomials) and not for
+standalone programs.
+
+AUTHOR:
+
+- Simon King, 2015-09-26
+
+diff --git a/src/cfinfo.c b/src/cfinfo.c
+index 293526b..9c1a004 100644
+--- a/src/cfinfo.c
++++ b/src/cfinfo.c
+@@ -215,7 +215,7 @@ int Lat_ReadInfo(Lat_Info *li, const char *basename)
+ 	    }
+ 	    for (i = 0; i < li->NCf; ++i)
+ 	    {
+-		ReadWord(f,&(li->Cf[i].idword),&(li->Cf[i].idpol),fn);
++		if (!ReadWord(f,&(li->Cf[i].idword),&(li->Cf[i].idpol),fn)) return -1;
+ 		if (StfMatch(f,i < li->NCf - 1 ? "," : "];") != 0)
+ 		{
+ 		    MTX_ERROR2("%s: %E",fn,MTX_ERR_FILEFMT);
+@@ -232,7 +232,7 @@ int Lat_ReadInfo(Lat_Info *li, const char *basename)
+ 	    }
+ 	    for (i = 0; i < li->NCf; ++i)
+ 	    {
+-		ReadWord(f,&(li->Cf[i].peakword),&(li->Cf[i].peakpol),fn);
++		if (!ReadWord(f,&(li->Cf[i].peakword),&(li->Cf[i].peakpol),fn)) return -1;
+ 		if (StfMatch(f,i < li->NCf - 1 ? "," : "];") != 0)
+ 		{
+ 		    MTX_ERROR2("%s: %E",fn,MTX_ERR_FILEFMT);
+diff --git a/src/chbasis.c b/src/chbasis.c
+index 34cf886..f1ee2e8 100644
+--- a/src/chbasis.c
++++ b/src/chbasis.c
+@@ -61,7 +61,8 @@ int MrChangeBasis(MatRep_t *rep, const Matrix_t *trans)
+ /** Conjugate a list @em gen of @em ngen square matrices over the same
+  *  field and of the same dimensions by a mattrix @em trans
+  *  and write the result into @em newgen. If @em gen == @em newgen, then
+- *  the previous content of @em newgen will be overridden. **/
++ *  the previous content of @em newgen will be overridden.
++ *  Return -1 on error and 0 on success. **/
+ int ChangeBasis(const Matrix_t *trans, int ngen, const Matrix_t *gen[],
+ 	Matrix_t *newgen[])
+ 
+@@ -83,18 +84,36 @@ int ChangeBasis(const Matrix_t *trans, int ngen, const Matrix_t *gen[],
+     }
+ 
+     Matrix_t *tmp = MatAlloc(trans->Field, trans->Nor, trans->Noc);
++    if (!tmp) return -1;
+     size_t tmpsize = FfCurrentRowSize*trans->Nor;
+     for (i = 0; i < ngen; ++i)
+     {
+         MTX_VERIFY(gen[i]->Nor==trans->Nor);
+         MTX_VERIFY(gen[i]->Noc==trans->Noc);
+         memset(tmp->Data, FF_ZERO, tmpsize);
+-        MatMulStrassen(tmp, trans, gen[i]);
++        if (!MatMulStrassen(tmp, trans, gen[i]))
++        {
++			MatFree(tmp);
++			return -1;
++		}
+         if ((const Matrix_t **)newgen == gen)
+             memset(newgen[i]->Data, FF_ZERO, tmpsize);
+         else
++        {
+             newgen[i] = MatAlloc(trans->Field, trans->Nor, trans->Noc);
+-        MatMulStrassen(newgen[i], tmp, bi);
++            if (!newgen[i])
++            {
++                MatFree(tmp);
++                MatFree(bi);
++                return -1;
++            }
++        }
++        if (!MatMulStrassen(newgen[i], tmp, bi))
++        {
++            MatFree(tmp);
++            MatFree(bi);
++            return -1;
++        }
+     }
+     MatFree(bi);
+     MatFree(tmp);
+diff --git a/src/ffio.c b/src/ffio.c
+index 92f9360..d2e3f1c 100644
+--- a/src/ffio.c
++++ b/src/ffio.c
+@@ -71,8 +71,11 @@ int FfReadRows(FILE *f, PTR buf, int n)
+         if (fread(b,FfTrueRowSize(FfNoc),1,f) != 1) break;
+ 	b += FfCurrentRowSize;
+     }
+-    if (ferror(f)) 
+-	MTX_ERROR("Read failed: %S");
++    if (ferror(f))
++    {
++        MTX_ERROR("Read failed: %S");
++        return -1;
++    }
+     return i;
+ }
+ 
+@@ -106,8 +109,11 @@ int FfWriteRows(FILE *f, PTR  buf, int n)
+         if (fwrite(b,FfTrueRowSize(FfNoc),1,f) != 1) break;
+ 	b += FfCurrentRowSize;
+     }
+-    if (ferror(f)) 
+-	MTX_ERROR("Write failed: %S");
++    if (ferror(f))
++    {
++        MTX_ERROR("Write failed: %S");
++        return -1;
++    }
+     return i;
+ }
+ 
+diff --git a/src/kernel-0.c b/src/kernel-0.c
+index 6ef2f72..431f01a 100644
+--- a/src/kernel-0.c
++++ b/src/kernel-0.c
+@@ -304,7 +304,10 @@ static FILE *OpenTableFile(int fl)
+     /* Create the table file.
+        ---------------------- */
+     if (FfMakeTables(fl) != 0)
+-	MTX_ERROR("Unable to build arithmetic tables");
++	{
++        MTX_ERROR("Unable to build arithmetic tables");
++        return NULL;
++    }
+     fd = SysFopen(fn,FM_READ|FM_LIB);
+     return fd;
+ }
+@@ -363,8 +366,7 @@ static int ReadTableFile(FILE *fd, int field)
+ 	return -1;
+     }
+     FfOrder = field;
+-    FfSetNoc(FfOrder);
+-    return 0;
++    return FfSetNoc(FfOrder);
+ }
+ 
+ 
+@@ -471,7 +473,7 @@ size_t FfTrueRowSize(int noc)
+  ** Embed a subfield.
+  ** @param a Element of the subfield field.
+  ** @param subfield Subfield order. Must be a divisor of the current field order.
+- ** @return @em a, embedded into the current field.
++ ** @return @em a, embedded into the current field, or 255 on error.
+  **/ 
+ 
+ FEL FfEmbed(FEL a, int subfield)
+@@ -482,7 +484,9 @@ FEL FfEmbed(FEL a, int subfield)
+ 	return a;
+     for (i = 0; mtx_embedord[i] != subfield && i < 4; ++i);
+     if (i >= 4)
+-	MTX_ERROR2("Cannot embed GF(%d) into GF(%d)",(int)subfield,(int)FfOrder);
++	{ MTX_ERROR2("Cannot embed GF(%d) into GF(%d)",(int)subfield,(int)FfOrder);
++      return (FEL)255;
++    }
+     return mtx_embed[i][a];
+ }
+ 
+@@ -498,6 +502,7 @@ FEL FfEmbed(FEL a, int subfield)
+  ** <tt>FfSetField(subfield)</tt>.
+  ** @param a Element of the current field.
+  ** @param subfield Subfield order. Must be a divisor of the current field order.
++ ** Return 255 on error.
+  **/
+ 
+ FEL FfRestrict(FEL a, int subfield)
+@@ -511,6 +516,7 @@ FEL FfRestrict(FEL a, int subfield)
+     {
+ 	MTX_ERROR2("Cannot restrict GF(%d) to GF(%d)",(int)FfOrder,
+ 	    (int)subfield);
++        return (FEL)255;
+     }
+     return mtx_restrict[i][a];
+ }
+diff --git a/src/maddmul.c b/src/maddmul.c
+index f5c171d..24ad3a5 100644
+--- a/src/maddmul.c
++++ b/src/maddmul.c
+@@ -59,7 +59,7 @@ Matrix_t *MatAddMul(Matrix_t *dest, const Matrix_t *src, FEL coeff)
+ 	   ------------ */
+ 	PTR dp = dest->Data, sp = src->Data;
+ 	int n;
+-	FfSetField(src->Field);
++	FfSetField(src->Field);  /* No error checking */
+ 	FfSetNoc(src->Noc);
+ 	for (n = src->Nor; n > 0; --n)
+ 	{
+diff --git a/src/maketabF.c b/src/maketabF.c
+index d7af83e..0fa26fb 100644
+--- a/src/maketabF.c
++++ b/src/maketabF.c
+@@ -175,7 +175,7 @@ static void polymod(POLY a, POLY b)
+    testprim() - Test for primitivity.
+    ----------------------------------------------------------------- */
+ 
+-static void testprim()
++static int testprim()
+ {
+     int i, a[256];
+ 
+@@ -187,7 +187,9 @@ static void testprim()
+ 	{
+ 	    fprintf(stderr,"*** a[%d]=%d.",i,a[i]);
+ 	    MTX_ERROR("Polynome is not primitive.");
++        return 1;
+ 	}
++    return 0;
+ }
+ 
+ 
+@@ -195,7 +197,7 @@ static void testprim()
+    initarith() - Initialize index and zech logarithm tables.
+    ----------------------------------------------------------------- */
+ 
+-static void initarith()
++static int initarith()
+ {	int i,elem;
+ 	POLY a;
+ 
+@@ -214,7 +216,7 @@ static void initarith()
+ 		polmultx(a);
+ 		polymod(a,irred);
+         }
+-	testprim();
++	if (testprim()) return 1;
+ 
+ 	/* Calculate zech logarithms
+ 	   ------------------------- */
+@@ -222,6 +224,7 @@ static void initarith()
+ 	{	elem = (int)((i%P)==P-1 ? i+1-P : i+1); /* add 1 */
+ 		zech[indx[i]]=indx[elem]; /* Zech-table=result */
+         }
++    return 0;
+ }
+ 
+ 
+@@ -314,7 +317,7 @@ static BYTE pack(BYTE a[8])
+ 	and initialize tables.
+    ----------------------------------------------------------------- */
+ 
+-static void writeheader()
++static int writeheader()
+ {
+     int i, j;
+ 
+@@ -324,6 +327,7 @@ static void writeheader()
+     {
+ 	perror(filename);
+ 	MTX_ERROR("Cannot open table file");
++    return 1;
+     }
+     for (CPM=1,maxmem=Q; (long)maxmem * Q <= 256L; ++CPM, maxmem *= Q);
+     for (i = 0; irrednrs[i] != (int) Q && irrednrs[i] != 0; ++i);
+@@ -333,7 +337,7 @@ static void writeheader()
+         for (j = 0; j <= MAXGRAD; j++)
+             irred[j] = irreducibles[i][MAXGRAD-j];
+ 	G = P;		/* Generator is X */
+-	initarith();	/* Init index- and Zech-tables */
++	if (initarith()) return 1;	/* Init index- and Zech-tables */
+     }
+     else
+     {	
+@@ -357,6 +361,7 @@ static void writeheader()
+     }
+     MESSAGE(1,("Generator   : %ld\n",info[1]));
+     MESSAGE(1,("Packing     : %ld/byte\n",info[3]));
++    return 0;
+ }
+ 
+ 
+@@ -364,14 +369,14 @@ static void writeheader()
+    checkq() - Set Q and N. Verify that Q is a prime power.
+    ----------------------------------------------------------------- */
+ 
+-static void checkq(long l)
++static int checkq(long l)
+ {
+     long q, d;
+ 
+     if (l < 2 || l > 256)
+     {
+-	fprintf(stderr,"Field order out of range (2-256)\n");
+-	exit(EXIT_ERR);
++	MTX_ERROR1("Field order out of range (2-256): %E", MTX_ERR_RANGE);
++	return 1;
+     }
+ 
+     Q = l;
+@@ -381,9 +386,10 @@ static void checkq(long l)
+        	q /= d;
+     if (q != 1)
+     {
+-	fprintf(stderr,"Illegal Field order\n");
+-	exit(EXIT_ERR);
++	MTX_ERROR("Illegal Field order\n");
++	return 1;
+     }
++    return 0;
+ }
+ 
+ 
+@@ -407,7 +413,7 @@ static void inittables()
+    mkembed() - Calculate embeddings of all subfields.
+    ----------------------------------------------------------------- */
+ 
+-static void mkembed()
++static int mkembed()
+ {
+     int n;	/* Degree of subfield over Z_p */
+     long q; /* subfield order */
+@@ -456,6 +462,7 @@ static void mkembed()
+ 	{
+ 	    fprintf(stderr,"*** q=%ld, Q=%ld.",q,Q);
+ 	    MTX_ERROR("Internal error.");
++        return 1;
+ 	}
+ 
+ 	/* Calculate a generator for the subfield
+@@ -502,13 +509,13 @@ static void mkembed()
+ 	    fflush(stdout);
+ 	}
+     }
++    return 0;
+ }
+ 
+ 
+ static int Init(int field)
+ {
+-    checkq(field);
+-    return 0;
++    return checkq(field);
+ }
+ 
+ /* -----------------------------------------------------------------
+@@ -526,7 +533,7 @@ int FfMakeTables(int field)
+        ---------- */
+     if (Init(field) != 0)
+ 	return 1;
+-    writeheader();			/* Open file and write header */
++    if (writeheader()) return 1;			/* Open file and write header */
+     inittables();
+ 
+     /* Make insert table
+@@ -618,7 +625,7 @@ int FfMakeTables(int field)
+ 	}
+     }
+ 
+-    mkembed();
++    if (mkembed()) return 1;
+ 
+     MESSAGE(1,("Writing tables to %s\n",filename));
+     if (
+@@ -639,6 +646,7 @@ int FfMakeTables(int field)
+     {
+ 	perror(filename);
+ 	MTX_ERROR("Error writing table file");
++    return 1;
+     }
+     fclose(fd);
+     return(0);
+diff --git a/src/matadd.c b/src/matadd.c
+index 54dbcb1..2d86d86 100644
+--- a/src/matadd.c
++++ b/src/matadd.c
+@@ -48,7 +48,7 @@ Matrix_t *MatAdd(Matrix_t *dest, const Matrix_t *src)
+        ------------------- */
+     dp = dest->Data;
+     sp = src->Data;
+-    FfSetField(src->Field);
++    FfSetField(src->Field);   /* No error checking */
+     FfSetNoc(src->Noc);
+     for (n = src->Nor; n > 0; --n)
+     {
+diff --git a/src/matclean.c b/src/matclean.c
+index e7307bf..16f02d6 100644
+--- a/src/matclean.c
++++ b/src/matclean.c
+@@ -53,7 +53,7 @@ int MatClean(Matrix_t *mat, const Matrix_t *sub)
+ 
+     /* Clean
+        ----- */
+-    FfSetNoc(mat->Noc);
++    FfSetNoc(mat->Noc);  /* No error checking */
+     for (i = 0; i < mat->Nor; ++i)
+     {
+ 	PTR m = MatGetPtr(mat,i);
+diff --git a/src/matcmp.c b/src/matcmp.c
+index b778ec4..d503285 100644
+--- a/src/matcmp.c
++++ b/src/matcmp.c
+@@ -38,7 +38,7 @@ MTX_DEFINE_FILE_INFO
+  ** not necessarily mean that an error has occured.
+  ** @param a First matrix.
+  ** @param b Second matrix.
+- ** @return 0 if the matrices are equal, nonzero otherwise (see description).
++ ** @return 0 if the matrices are equal, nonzero otherwise (see description), -2 on error.
+  **/
+ 
+ int MatCompare(const Matrix_t *a, const Matrix_t *b)
+@@ -50,7 +50,7 @@ int MatCompare(const Matrix_t *a, const Matrix_t *b)
+     if (!MatIsValid(a) || !MatIsValid(b))
+     {
+ 	MTX_ERROR1("%E",MTX_ERR_BADARG);
+-	return -1;
++	return -2;
+     }
+ 
+     /* Compare fields and dimensions
+@@ -65,7 +65,7 @@ int MatCompare(const Matrix_t *a, const Matrix_t *b)
+     /* Compare the entries row by row. We do not use memcmp on the
+        whole matrix because we must ignore padding bytes.
+        ----------------------------------------------------------- */
+-    FfSetField(a->Field);
++    FfSetField(a->Field);  /* No error checking */
+     FfSetNoc(a->Noc);
+     for (i = 0; i < a->Nor; ++i)
+     {
+diff --git a/src/matcopy.c b/src/matcopy.c
+index 75852dd..c3e7850 100644
+--- a/src/matcopy.c
++++ b/src/matcopy.c
+@@ -105,7 +105,7 @@ int MatCopyRegion(Matrix_t *dest, int destrow, int destcol,
+ 	{
+ #ifdef PARANOID
+ 	    FEL f;
+-	    FfSetNoc(src->Noc);
++	    FfSetNoc(src->Noc);  /* No error checking */
+ 	    f = FfExtract(s,k);
+ 	    FfSetNoc(dest->Noc);
+ 	    FfInsert(d,destcol+k-col1,f);
+diff --git a/src/matcore.c b/src/matcore.c
+index 1f27dfd..0dc9d92 100644
+--- a/src/matcore.c
++++ b/src/matcore.c
+@@ -131,7 +131,7 @@ Matrix_t *MatAlloc(int field, int nor, int noc)
+ 	SysFree(m);
+ 	return NULL;
+     }
+-    FfSetNoc(noc);
++    if (FfSetNoc(noc)) return NULL;
+     m->Magic = MAT_MAGIC;
+     m->Field = field;
+     m->Nor = nor;
+diff --git a/src/matcut.c b/src/matcut.c
+index fde0662..f274311 100644
+--- a/src/matcut.c
++++ b/src/matcut.c
+@@ -79,11 +79,12 @@ Matrix_t *MatCut(const Matrix_t *src, int row1, int col1, int nrows, int ncols)
+     /* Initialize pointers to the source and destination matrix 
+        -------------------------------------------------------- */
+     s = MatGetPtr(src,row1);
++    if (!s) return NULL;
+     d = result->Data;
+ 
+     /* Copy the requested data
+        ----------------------- */
+-    FfSetNoc(ncols);
++    if (FfSetNoc(ncols)) return NULL;
+     for (n = nrows; n > 0; --n)
+     {
+ 	if (col1 == 0)
+@@ -95,9 +96,9 @@ Matrix_t *MatCut(const Matrix_t *src, int row1, int col1, int nrows, int ncols)
+ 	    {
+ #ifdef PARANOID
+ 		FEL f;
+-		FfSetNoc(src->Noc);
++		FfSetNoc(src->Noc);  /* No error checking */
+ 		f = FfExtract(s,col1+k);
+-		FfSetNoc(ncols);
++		FfSetNoc(ncols);  /* error was checked above */
+ 		FfInsert(d,k,f);
+ #else
+ 		FfInsert(d,k,FfExtract(s,col1+k));
+diff --git a/src/matech.c b/src/matech.c
+index ed31cf4..ee52ebe 100644
+--- a/src/matech.c
++++ b/src/matech.c
+@@ -124,7 +124,7 @@ int MatEchelonize(Matrix_t *mat)
+ 
+     /* Build the pivot table
+        --------------------- */
+-    FfSetField(mat->Field);
++    FfSetField(mat->Field);  /* No error checking */
+     FfSetNoc(mat->Noc);
+     rank = zmkechelon(mat->Data,mat->Nor,mat->Noc,mat->PivotTable,is_pivot);
+ 
+@@ -163,13 +163,14 @@ long MatNullity(const Matrix_t *mat)
+  ** This function calculates the dimension of the null-space of a matrix
+  ** and deletes the matrix.
+  ** @param mat Pointer to the matrix.
+- ** @return Nullity of @em mat, or -$ on error.
++ ** @return Nullity of @em mat, or $-1$ on error.
+  **/
+ 
+ long MatNullity__(Matrix_t *mat)
+ {
+     long nul;
+-    MatEchelonize(mat);
++    if (!mat) return -1;
++    if (MatEchelonize(mat)==-1) return -1;
+     nul = mat->Noc - mat->Nor;
+     MatFree(mat);
+     return nul;
+diff --git a/src/matins.c b/src/matins.c
+index 45c31e4..50fe9c1 100644
+--- a/src/matins.c
++++ b/src/matins.c
+@@ -54,7 +54,7 @@ Matrix_t *MatInsert_(Matrix_t *mat, const Poly_t *pol)
+ 	return NULL;
+     }
+ 
+-    FfSetField(mat->Field);
++    FfSetField(mat->Field);  /* No error checking */
+     FfSetNoc(nor);
+ 
+     /* Special case: p(x) = 0
+@@ -81,7 +81,10 @@ Matrix_t *MatInsert_(Matrix_t *mat, const Poly_t *pol)
+     /* Evaluate p(A)
+        ------------- */
+     if (pol->Degree > 1) 
+-	x = MatDup(mat);
++	{
++		x = MatDup(mat);
++		if (!x) return NULL;
++	}
+     if ((f = pol->Data[pol->Degree]) != FF_ONE)
+     {
+ 	for (l = nor, v = mat->Data; l > 0; --l, FfStepPtr(&v))
+@@ -147,6 +150,7 @@ Matrix_t *MatInsert(const Matrix_t *mat, const Poly_t *pol)
+     if (pol->Degree == 0)
+     {
+ 	x = MatAlloc(mat->Field,nor,nor);
++    if (!x) return NULL;
+ 	for (l = 0, v = x->Data; l < nor; ++l, FfStepPtr(&v))
+ 	    FfInsert(v,l,pol->Data[0]);
+ 	return x;
+@@ -155,6 +159,7 @@ Matrix_t *MatInsert(const Matrix_t *mat, const Poly_t *pol)
+     /* Evaluate p(A)
+        ------------- */
+     x = MatDup(mat);
++    if (!x) return NULL;
+     if ((f = pol->Data[pol->Degree]) != FF_ONE)
+     {
+ 	for (l = nor, v = x->Data; l > 0; --l, FfStepPtr(&v))
+diff --git a/src/matinv.c b/src/matinv.c
+index 217eb0e..990cbe7 100644
+--- a/src/matinv.c
++++ b/src/matinv.c
+@@ -114,6 +114,7 @@ Matrix_t *MatInverse(const Matrix_t *mat)
+     /* Copy matrix into workspace
+        -------------------------- */
+     tmp = FfAlloc(mat->Nor);
++    if (!tmp) return NULL;
+     memcpy(tmp,mat->Data,FfCurrentRowSize * mat->Nor);
+ 
+     /* Inversion
+diff --git a/src/matmul.c b/src/matmul.c
+index 20f5e88..bed30fc 100644
+--- a/src/matmul.c
++++ b/src/matmul.c
+@@ -63,7 +63,7 @@ Matrix_t *MatMul(Matrix_t *dest, const Matrix_t *src)
+ 
+     /* Matrix multiplication
+        --------------------- */
+-    FfSetField(src->Field);
++    FfSetField(src->Field); /* no error checking, since the matrix *exists* */
+     FfSetNoc(src->Noc);
+     result = tmp = FfAlloc(dest->Nor);
+     if (result == NULL)
+diff --git a/src/matnull.c b/src/matnull.c
+index 4f28566..2550b96 100644
+--- a/src/matnull.c
++++ b/src/matnull.c
+@@ -27,6 +27,8 @@ MTX_DEFINE_FILE_INFO
+     - |piv| contains a pivot table for the null space.
+     If |flags| is nonzero, the null-space is not reduced to echelon form,
+     and the contents of |piv| are undefined.
++
++    Return -1 on error, the dimension of the null-space on success.
+  ** @see 
+  **/
+ 
+@@ -40,7 +42,7 @@ static long znullsp(PTR matrix, long nor, int *piv, PTR nsp, int flags)
+ 
+     /* Make the identity matrix in <nsp>.
+        ---------------------------------- */
+-    FfSetNoc(nor);
++    if (FfSetNoc(nor)) return -1;
+     x = nsp;
+     for (i = 0; i < nor; ++i)
+     {
+@@ -61,13 +63,12 @@ static long znullsp(PTR matrix, long nor, int *piv, PTR nsp, int flags)
+ 
+ 	for (k = 0; k < i; ++k)
+ 	{
+-	    FfSetNoc(noc);
++	    FfSetNoc(noc);  /* No error checking, since noc used to be the previously assigned number of columns */
+ 	    if ((p = piv[k]) >= 0 && (f = FfExtract(x,p)) != FF_ZERO)
+ 	    {
+ 		f = FfNeg(FfDiv(f,FfExtract(xx,p)));
+-		FfSetNoc(noc);
+ 		FfAddMulRow(x,xx,f);
+-		FfSetNoc(nor);
++		FfSetNoc(nor);  /* we have asserted above that it doesn't fail */
+ 		FfAddMulRow(y,yy,f);
+ 	    }
+ 	    FfSetNoc(noc);
+@@ -151,11 +152,21 @@ Matrix_t *MatNullSpace_(Matrix_t *mat, int flags)
+     if (nsp == NULL) 
+ 	return NULL;
+     nsp->PivotTable = NREALLOC(nsp->PivotTable,int,mat->Nor);
++    if (!nsp->PivotTable)
++    {
++        MatFree(nsp);
++        return NULL;
++    }
+ 
+     /* Calculate the null-space
+        ------------------------ */
+-    FfSetNoc(mat->Noc);
++    FfSetNoc(mat->Noc);  /* No error checking */
+     dim = znullsp(mat->Data,mat->Nor,nsp->PivotTable,nsp->Data,flags);
++    if (dim==-1)
++    {
++        MatFree(nsp);
++        return NULL;
++    }
+     if (flags)
+     {
+ 	SysFree(nsp->PivotTable);
+diff --git a/src/matorder.c b/src/matorder.c
+index 16aec74..24b31a3 100644
+--- a/src/matorder.c
++++ b/src/matorder.c
+@@ -32,7 +32,7 @@ MTX_DEFINE_FILE_INFO
+  ** the order is greater than 1000000, or if the order on any cyclic
+  ** subspace is greater than 1000.
+  ** @param mat Pointer to the matrix.
+- ** @return The order of @em mat, or 1 on error.
++ ** @return The order of @em mat, or -1 on error.
+  **/
+ 
+ int MatOrder(const Matrix_t *mat)
+@@ -59,15 +59,29 @@ int MatOrder(const Matrix_t *mat)
+     FfSetNoc(mat->Noc);
+     nor = mat->Nor;
+     m1 = FfAlloc(nor);
++    if (!m1) return -1;
+     memcpy(m1,mat->Data,FfCurrentRowSize * nor);
+     bend = basis = FfAlloc(nor+1);
++    if (!bend)
++    {
++        SysFree(m1);
++        return -1;
++    }
+ 
+     piv = NALLOC(int,nor+1);
+     done = NALLOC(char,nor);
++    if (!piv || !done)
++    { SysFree(m1);
++      return -1;
++    }
+     memset(done,0,(size_t)nor);
+     v1 = FfAlloc(1);
+     v2 = FfAlloc(1);
+     v3 = FfAlloc(1);
++    if (!v1 || !v2 || !v3)
++    { SysFree(m1);
++      return -1;
++    }
+     tord = ord = 1;
+     dim = 0;
+     j1 = 1;
+diff --git a/src/matpivot.c b/src/matpivot.c
+index abe342a..c843282 100644
+--- a/src/matpivot.c
++++ b/src/matpivot.c
+@@ -71,7 +71,6 @@ static int zmkpivot(PTR matrix, int nor, int noc, int *piv, int *ispiv)
+ 
+ int MatPivotize(Matrix_t *mat)
+ {
+-    int rc;
+     int *newtab;
+     static int *is_pivot = NULL;
+     static int maxnoc = -1;
+@@ -106,9 +105,7 @@ int MatPivotize(Matrix_t *mat)
+        --------------------- */
+     FfSetField(mat->Field);
+     FfSetNoc(mat->Noc);
+-    rc = zmkpivot(mat->Data,mat->Nor,mat->Noc,mat->PivotTable,is_pivot);
+-
+-    return rc;
++    return zmkpivot(mat->Data,mat->Nor,mat->Noc,mat->PivotTable,is_pivot);
+ }
+ 
+ /**
+diff --git a/src/matpwr.c b/src/matpwr.c
+index cd66b5e..b06e5b2 100644
+--- a/src/matpwr.c
++++ b/src/matpwr.c
+@@ -119,8 +119,14 @@ Matrix_t *MatPower(const Matrix_t *mat, long n)
+     FfSetField(mat->Field);
+     FfSetNoc(mat->Noc);
+     tmp = FfAlloc(FfNoc);
++    if (!tmp) return NULL;
+     memcpy(tmp,mat->Data,FfCurrentRowSize * FfNoc);
+     tmp2 = FfAlloc(FfNoc);
++    if (!tmp2)
++    {
++        SysFree(tmp);
++        return NULL;
++    }
+     result = MatAlloc(mat->Field,mat->Nor,mat->Noc);
+     if (result != NULL)
+ 	matpwr_(n,tmp,result->Data,tmp2);
+diff --git a/src/matread.c b/src/matread.c
+index 031d100..06e6e6b 100644
+--- a/src/matread.c
++++ b/src/matread.c
+@@ -46,8 +46,9 @@ Matrix_t *MatRead(FILE *f)
+ 	return NULL;
+     if (FfReadRows(f,m->Data,m->Nor) != m->Nor)
+     {
+-	MatFree(m);
+-	return NULL;
++        MTX_ERROR("Number of given rows does not coincide with given row number");
++        MatFree(m);
++        return NULL;
+     }
+     return m;
+ }
+diff --git a/src/mattrace.c b/src/mattrace.c
+index f500248..772a6e4 100644
+--- a/src/mattrace.c
++++ b/src/mattrace.c
+@@ -21,7 +21,7 @@
+  ** This function calculates the sum of all diagonal elements of a matrix.
+  ** Note that the matrix need not be square.
+  ** @param mat Pointer to the matrix.
+- ** @return Trace of @a mat, @c FF_ZERO on error.
++ ** @return Trace of @a mat, @c 255 on error.
+  **/
+ 
+ FEL MatTrace(const Matrix_t *mat)
+@@ -35,7 +35,7 @@ FEL MatTrace(const Matrix_t *mat)
+        ------------------ */
+ #ifdef DEBUG
+     if (!MatIsValid(mat))
+-	return FF_ZERO;
++	return (FEL)255;
+ #endif
+ 
+     maxi = mat->Nor > mat->Noc ? mat->Noc : mat->Nor;
+diff --git a/src/matwrite.c b/src/matwrite.c
+index 1fb6af3..b364e80 100644
+--- a/src/matwrite.c
++++ b/src/matwrite.c
+@@ -44,7 +44,10 @@ int MatWrite(const Matrix_t *mat, FILE *f)
+     FfSetField(mat->Field);
+     FfSetNoc(mat->Noc);
+     if (FfWriteRows(f,mat->Data,mat->Nor) != mat->Nor)
+-	return -1;
++	{
++        MTX_ERROR("Cannot write rows");
++        return -1;
++    }
+     return 0;
+ }
+ 
+@@ -75,7 +78,10 @@ int MatSave(const Matrix_t *mat, const char *fn)
+     i = MatWrite(mat,f);
+     fclose(f);
+     if (i != 0)
+-	MTX_ERROR1("Cannot write matrix to %s",fn);
++	{
++        MTX_ERROR1("Cannot write matrix to %s",fn);
++        return -1;
++    }
+     return i;
+ }
+ 
+diff --git a/src/meataxe.h b/src/meataxe.h
+index 368b37b..0efa7dd 100644
+--- a/src/meataxe.h
++++ b/src/meataxe.h
+@@ -135,8 +135,8 @@ PTR FfSubRowPartialReverse(PTR dest, PTR src, int first, int len);
+ PTR FfAlloc(int nor);
+ int FfCmpRows(PTR p1, PTR p2);
+ void FfCleanRow(PTR row, PTR matrix, int nor, const int *piv);
+-void FfCleanRow2(PTR row, PTR matrix, int nor, const int *piv, PTR row2);
+-void FfCleanRowAndRepeat(PTR row, PTR mat, int nor, const int *piv, 
++int FfCleanRow2(PTR row, PTR matrix, int nor, const int *piv, PTR row2);
++int FfCleanRowAndRepeat(PTR row, PTR mat, int nor, const int *piv, 
+     PTR row2, PTR mat2);
+ void FfCopyRow(PTR dest, PTR src);
+ FEL FfEmbed(FEL a, int subfield);
+diff --git a/src/mmulscal.c b/src/mmulscal.c
+index 281be16..9bff3ff 100644
+--- a/src/mmulscal.c
++++ b/src/mmulscal.c
+@@ -21,7 +21,7 @@
+  ** Multiply a Matrix by a Constant.
+  ** @param dest Pointer to the matrix.
+  ** @param coeff Value to multiply with.
+- ** @return The function returns @a dest.
++ ** @return The function returns @a dest, or NULL on error in debug mode only.
+  **/
+ 
+ Matrix_t *MatMulScalar(Matrix_t *dest, FEL coeff)
+diff --git a/src/mtensor.c b/src/mtensor.c
+index b2b9ac7..2fb90ca 100644
+--- a/src/mtensor.c
++++ b/src/mtensor.c
+@@ -85,6 +85,11 @@ Matrix_t *MatTensor(const Matrix_t *m1, const Matrix_t *m2)
+ 	   ---------------------------------------- */
+ 	x1 = m1->Data;
+ 	x3 = MatGetPtr(temat,i2);
++    if (!x3)
++    {
++        MatFree(temat);
++        return NULL;
++    }
+ 	FfSetNoc(temat->Noc);
+ 
+ 	/* Loop through all rows of <m1>
+diff --git a/src/quotient.c b/src/quotient.c
+index f44b556..eea0754 100644
+--- a/src/quotient.c
++++ b/src/quotient.c
+@@ -82,16 +82,23 @@ Matrix_t *QProjection(const Matrix_t *subspace, const Matrix_t *vectors)
+     sdim = subspace->Nor;
+     qdim = subspace->Noc - sdim;
+     result = MatAlloc(subspace->Field,vectors->Nor,qdim);
++    if (!result) return NULL;
+ 
+     /* Calculate the projection
+        ------------------------ */
+     FfSetNoc(subspace->Noc);
+     tmp = FfAlloc(1);
++    if (!tmp) return NULL;
+     non_piv = subspace->PivotTable + subspace->Nor;
+     for (i = 0; i < vectors->Nor; ++i)
+     {
+ 	int k;
+ 	PTR q = MatGetPtr(result,i);
++    if (!q)
++    {
++        SysFree(tmp);
++        return NULL;
++    }
+ 	FfCopyRow(tmp,MatGetPtr(vectors,i));
+ 	FfCleanRow(tmp,subspace->Data,sdim,subspace->PivotTable);
+ 	for (k = 0; k < qdim; ++k)
+@@ -158,14 +165,20 @@ Matrix_t *QAction(const Matrix_t *subspace, const Matrix_t *gen)
+ 
+     /* Calculate the action on the quotient
+        ------------------------------------ */
+-    FfSetNoc(dim);
++    FfSetNoc(dim);  /* No error checking, since dim is the ->Noc of an existing matrix */
+     tmp = FfAlloc(1);
++    if (!tmp) return NULL;
+     piv = subspace->PivotTable;
+     non_piv = piv + subspace->Nor;
+     for (k = 0; k < qdim; ++k)
+     {
+ 	int l;
+ 	PTR qx = MatGetPtr(action,k);
++    if (!qx)
++    {
++        SysFree(tmp);
++        return NULL;
++    }
+ 	FfCopyRow(tmp,MatGetPtr(gen,non_piv[k]));
+ 	FfCleanRow(tmp,subspace->Data,sdim,piv);
+ 	for (l = 0; l < qdim; ++l)
+diff --git a/src/saction.c b/src/saction.c
+index adae3cf..0aba44d 100644
+--- a/src/saction.c
++++ b/src/saction.c
+@@ -68,8 +68,14 @@ Matrix_t *SAction(const Matrix_t *subspace, const Matrix_t *gen)
+     sdim = subspace->Nor;
+     FfSetField(subspace->Field);
+     action = MatAlloc(FfOrder,sdim,sdim);
+-    FfSetNoc(dim);
++    if (!action) return NULL;
++    FfSetNoc(dim);  /* No error checking, since dim is the ->Noc of an existing matrix */
+     tmp = FfAlloc(1);
++    if (!tmp)
++    {
++        MatFree(action);
++        return NULL;
++    }
+ 
+     /* Calaculate the action.
+        ---------------------- */
+@@ -77,6 +83,12 @@ Matrix_t *SAction(const Matrix_t *subspace, const Matrix_t *gen)
+     {
+ 	PTR xi = MatGetPtr(subspace,i);
+ 	PTR yi = MatGetPtr(action,i);
++    if (!xi || !yi)
++    {
++        MatFree(action);
++        SysFree(tmp);
++        return NULL;
++    }
+ 	FEL f;
+ 
+ 	/* Calculate the image of the <i>-th row of <subspace>.
+@@ -85,10 +97,20 @@ Matrix_t *SAction(const Matrix_t *subspace, const Matrix_t *gen)
+ 
+ 	/* Clean the image with the subspace and store coefficients.
+ 	   --------------------------------------------------------- */
+-	FfCleanRow2(tmp,subspace->Data,sdim,subspace->PivotTable,yi);
++	if (FfCleanRow2(tmp,subspace->Data,sdim,subspace->PivotTable,yi))
++    {
++        MatFree(action);
++        SysFree(tmp);
++        return NULL;
++    }
+ 	if (FfFindPivot(tmp,&f) >= 0)
+-	    MTX_ERROR("Split(): Subspace not invariant");
++	{
++        MatFree(action);
++        SysFree(tmp);
++        MTX_ERROR("Split(): Subspace not invariant");
++        return NULL;
+     }
++	}
+ 
+     /* Clean up and return the result.
+        ------------------------------- */
+diff --git a/src/stabpwr.c b/src/stabpwr.c
+index ff33bc6..01282ab 100644
+--- a/src/stabpwr.c
++++ b/src/stabpwr.c
+@@ -68,15 +68,18 @@ int StablePower_(Matrix_t *mat, int *pwr, Matrix_t **ker)
+        --------------------------- */
+     p = 1;
+     k1 = MatNullSpace(mat);
+-    MatMul(mat,mat);
++    if (!k1) return -1;
++    if (!MatMul(mat,mat)) return -1;
+     k2 = MatNullSpace(mat);
++    if (!k2) return -1;
+     while (k2->Nor > k1->Nor)
+     {
+ 	p *= 2;
+ 	MatFree(k1);
+ 	k1 = k2;
+-	MatMul(mat,mat);
++	if (!MatMul(mat,mat)) return -1;
+ 	k2 = MatNullSpace(mat);
++    if (!k2) return -1;
+     }
+     MatFree(k2);
+ 
+diff --git a/src/sumint.c b/src/sumint.c
+index 278acd8..905fa79 100644
+--- a/src/sumint.c
++++ b/src/sumint.c
+@@ -77,7 +77,7 @@ int FfSumAndIntersection(PTR wrk1, int *nor1, int *nor2, PTR wrk2, int *piv)
+     {
+ 	FEL f;
+ 	int p;
+-	FfCleanRowAndRepeat(x1,wrk1,k,piv,x2,wrk2);
++	if (FfCleanRowAndRepeat(x1,wrk1,k,piv,x2,wrk2)) return -1;
+ 	if ((p = FfFindPivot(x1,&f)) < 0)
+ 	    continue;	/* Null row - ignore */
+ 	if (k < i)
+diff --git a/src/temap.c b/src/temap.c
+index 7ba445a..4c2d493 100644
+--- a/src/temap.c
++++ b/src/temap.c
+@@ -74,17 +74,21 @@ Matrix_t *TensorMap(Matrix_t *vec, const Matrix_t *a, const Matrix_t *b)
+     for (i = 0; i < vec->Nor; ++i)
+     {
+ 	Matrix_t *tmp = MatTransposed(a);
++    if (!tmp) return NULL;
+ 	Matrix_t *v = VectorToMatrix(vec,i,b->Nor);
+ 	if (v == NULL)
+ 	{
+ 	    MTX_ERROR("Conversion failed");
+-	    break;
++	    return NULL;
+ 	}
+-	MatMul(tmp,v);
++	if (!MatMul(tmp,v)) return NULL;
+ 	MatFree(v);
+-	MatMul(tmp,b);
++	if (!MatMul(tmp,b)) return NULL;
+ 	if (MatrixToVector(tmp,result,i))
+-	    MTX_ERROR("Conversion failed");
++    {
++        MTX_ERROR("Conversion failed");
++        return NULL;
++    }
+ 	MatFree(tmp);
+     }
+     return result;
+diff --git a/src/vec2mat.c b/src/vec2mat.c
+index 1047805..e76ad88 100644
+--- a/src/vec2mat.c
++++ b/src/vec2mat.c
+@@ -63,8 +63,11 @@ Matrix_t *VectorToMatrix(Matrix_t *vecs, int n, int noc)
+ 	return NULL;
+     for (i = 0; i < result->Nor; ++i)
+     {
+-	if (MatCopyRegion(result,i,0, vecs,n,i*noc,1,noc) != 0)
+-	    MTX_ERROR("Copy failed");
++        if (MatCopyRegion(result,i,0, vecs,n,i*noc,1,noc) != 0)
++        {
++            MTX_ERROR("Copy failed");
++            return NULL;
++        }
+     }
+     return result;
+ }
+diff --git a/src/window.c b/src/window.c
+index 9c87694..fbeb943 100644
+--- a/src/window.c
++++ b/src/window.c
+@@ -69,7 +69,11 @@ MatrixWindow_t *WindowAlloc(int fl, int nor, size_t rowsize)
+         MTX_ERROR1("%E",MTX_ERR_NOMEM);
+         return NULL;
+     }
+-    FfSetField(fl);
++    if (FfSetField(fl))
++    {
++        free(out);
++        return NULL;
++    }
+     out->Matrix = MatAlloc(fl, nor, rowsize*sizeof(long)*MPB);
+     if (out->Matrix == NULL)
+     {
+@@ -266,7 +270,8 @@ __asm__("    popl %ebx\n"
+ 
+ /** dest := left+right
+    left and right must be distinct, but one of them may coincide with dest -- under the assumption
+-   that, in that case, the ambient matrices coincide as well. **/
++   that, in that case, the ambient matrices coincide as well.
++   Return dest, or NULL on error (the only error may occur in a compatibility check). **/
+ MatrixWindow_t *WindowSum(MatrixWindow_t *dest, MatrixWindow_t *left, MatrixWindow_t *right)
+ {
+   PTR x, result, tmp;
+@@ -335,6 +340,7 @@ MatrixWindow_t *WindowSum(MatrixWindow_t *dest, MatrixWindow_t *left, MatrixWind
+ /** dest := left-right
+    left and right must be distinct, but one of them may coincide with dest -- under the assumption
+    that, in that case, the ambient matrices coincide as well.
++   Return dest, or NULL on error (the only error may occur in a compatibility check).
+ **/
+ MatrixWindow_t *WindowDif(MatrixWindow_t *dest, MatrixWindow_t *left, MatrixWindow_t *right)
+ {
+@@ -407,7 +413,7 @@ MatrixWindow_t *WindowDif(MatrixWindow_t *dest, MatrixWindow_t *left, MatrixWind
+    can write the result into it. Moreover, the chunk of memory pointed at by dest MUST be disjoint
+    from the chunks for left and right!
+ 
+-   Dimensions are not tested!
++   Dimensions are not tested, always dest will be returned (no error value).
+ **/
+ MatrixWindow_t *WindowAddMul(MatrixWindow_t *dest, MatrixWindow_t *left, MatrixWindow_t *right)
+ {
+@@ -617,7 +623,7 @@ int StrassenStep(MatrixWindow_t *dest_win, MatrixWindow_t *A_win, MatrixWindow_t
+   S2->RowSize = A_sub_rowsize;
+   S2->Matrix = X->Matrix;
+   S2->ULCorner = X->ULCorner;
+-  WindowDif(S2, A00, A10);
++  WindowDif(S2, A00, A10); /* No error checking, as we know that the windows are compatible */
+   /*
+   printf("1.  S2 = A00-A10 in X\n");
+   WindowShow(X);
+@@ -653,7 +659,7 @@ int StrassenStep(MatrixWindow_t *dest_win, MatrixWindow_t *A_win, MatrixWindow_t
+   S0->RowSize = A_sub_rowsize;
+   S0->Matrix = X->Matrix;
+   S0->ULCorner = X->ULCorner;
+-  WindowSum(S0, A10, A11);
++  WindowSum(S0, A10, A11); /* no error checking here and below, as we know the dimensions of the windows */
+   /*
+   printf("4.  S0 = A10+A11 in X\n");
+   WindowShow(X);
+diff --git a/src/zcleanrow.c b/src/zcleanrow.c
+index b4dcb30..d36a165 100644
+--- a/src/zcleanrow.c
++++ b/src/zcleanrow.c
+@@ -63,10 +63,10 @@ void FfCleanRow(PTR row, PTR matrix, int nor, const int *piv)
+  ** @param nor Number of rows.
+  ** @param piv Pivot table for @em matrix.
+  ** @param row2 Pointer to row where the operations are recorded.
+- ** @return Always 0.
++ ** @return 0, or 1 on error.
+  **/
+ 
+-void FfCleanRow2(PTR row, PTR mat, int nor, const int *piv, PTR row2)
++int FfCleanRow2(PTR row, PTR mat, int nor, const int *piv, PTR row2)
+ {
+     int i;
+     PTR x;
+@@ -74,7 +74,7 @@ void FfCleanRow2(PTR row, PTR mat, int nor, const int *piv, PTR row2)
+     if (row2 == NULL || piv == NULL)
+     {
+ 	MTX_ERROR1("%E",MTX_ERR_BADARG);
+-	return;
++	return 1;
+     }
+     for (i = 0, x = mat; i < nor; ++i, FfStepPtr(&x))
+     {
+@@ -86,6 +86,7 @@ void FfCleanRow2(PTR row, PTR mat, int nor, const int *piv, PTR row2)
+ 	    FfInsert(row2,i,f);
+ 	}
+     }
++    return 0;
+ }
+ 
+ 
+@@ -100,10 +101,10 @@ void FfCleanRow2(PTR row, PTR mat, int nor, const int *piv, PTR row2)
+  ** @param piv Pivot table for @em mat.
+  ** @param row2 Pointer to the second row to be cleaned.
+  ** @param mat2 Matrix to the second matrix.
+- ** @return Always 0.
++ ** @return 0, or 1 on error.
+  **/
+ 
+-void FfCleanRowAndRepeat(PTR row, PTR mat, int nor, const int *piv, PTR row2, PTR mat2)
++int FfCleanRowAndRepeat(PTR row, PTR mat, int nor, const int *piv, PTR row2, PTR mat2)
+ {
+     int i;
+     PTR x, x2;
+@@ -112,7 +113,7 @@ void FfCleanRowAndRepeat(PTR row, PTR mat, int nor, const int *piv, PTR row2, PT
+     if (row2 == NULL || piv == NULL || row2 == NULL || mat2 == NULL)
+     {
+ 	MTX_ERROR1("%E",MTX_ERR_BADARG);
+-	return;
++	return 1;
+     }
+ #endif
+     for (i = 0, x = mat, x2 = mat2; i < nor; ++i, FfStepPtr(&x), FfStepPtr(&x2))
+@@ -125,6 +126,7 @@ void FfCleanRowAndRepeat(PTR row, PTR mat, int nor, const int *piv, PTR row2, PT
+ 	    FfAddMulRow(row2,x2,f);
+ 	}
+     }
++    return 0;
+ }
+ 
+ 
diff --git a/src/sage/libs/meataxe.pxd b/src/sage/libs/meataxe.pxd
index fc76bfc781e..68878f3fa19 100644
--- a/src/sage/libs/meataxe.pxd
+++ b/src/sage/libs/meataxe.pxd
@@ -56,8 +56,8 @@ cdef extern from "meataxe.h":
     # FEL FfMul(FEL a, FEL b)
     # FEL FfDiv(FEL a, FEL b)
     # FEL FfInv(FEL a)
-    # FEL FfEmbed(FEL a, int subfield)
-    # FEL FfRestrict(FEL a, int subfield)
+    # FEL FfEmbed(FEL a, int subfield) except 255
+    # FEL FfRestrict(FEL a, int subfield) except 255
     FEL FfFromInt(int l)
     int FfToInt(FEL f)
 
@@ -97,8 +97,8 @@ cdef extern from "meataxe.h":
     ## Basic memory operations
     Matrix_t *MatAlloc(int field, int nor, int noc) except NULL
     int MatFree(Matrix_t *mat)
-    PTR MatGetPtr(Matrix_t *mat, int row)
-    int MatCompare(Matrix_t *a, Matrix_t *b) except? -1
+    PTR MatGetPtr(Matrix_t *mat, int row) except NULL
+    int MatCompare(Matrix_t *a, Matrix_t *b) except -2
     int MatCopyRegion(Matrix_t *dest, int destrow, int destcol, Matrix_t *src, int row1, int col1, int nrows, int ncols) except -1
     Matrix_t *MatCut(Matrix_t *src, int row1, int col1, int nrows, int ncols) except NULL
     Matrix_t *MatCutRows(Matrix_t *src, int row1, int nrows) except NULL
@@ -115,11 +115,15 @@ cdef extern from "meataxe.h":
     Matrix_t *MatMul(Matrix_t *dest, Matrix_t *src) except NULL
     Matrix_t *MatMulScalar(Matrix_t *dest, FEL coeff) except NULL
     Matrix_t *MatPower(Matrix_t *mat, long n) except NULL
-    FEL MatTrace(Matrix_t *mat)
+    int StablePower(Matrix_t *mat, int *pwr, Matrix_t **ker) except -1
+    FEL MatTrace(Matrix_t *mat) except 255
     Matrix_t *MatMulStrassen(Matrix_t *dest, Matrix_t *A, Matrix_t *B) except NULL
     void StrassenSetCutoff(size_t size)
 
     ## "Higher" Arithmetic
+    Matrix_t *MatTensor(Matrix_t *m1, Matrix_t *m2) except NULL
+    Matrix_t *TensorMap(Matrix_t *vec, Matrix_t *a, Matrix_t *b) except NULL
+    
     int MatClean(Matrix_t *mat, Matrix_t *sub) except -1
     int MatEchelonize(Matrix_t *mat) except -1
     int MatOrder(Matrix_t *mat) except? -1