diff --git a/library/matrix/matrix.c b/library/matrix/matrix.c
index 5c06bdc8a..0ee9f68af 100644
--- a/library/matrix/matrix.c
+++ b/library/matrix/matrix.c
@@ -1,22 +1,20 @@
 /*************************************************************************
-*									 *
-*	 YAP Prolog 							 *
-*									 *
-*	Yap Prolog was developed at NCCUP - Universidade do Porto	 *
-*									 *
-* Copyright L.Damas, V.S.Costa and Universidade do Porto 1985-1997	 *
-*									 *
-**************************************************************************
-*									 *
-* File:		matrix.c						 *
-* Last rev:								 *
-* mods:									 *
-* comments:	numerical arrays		                         *
-*									 *
-*************************************************************************/
+ *									 *
+ *	 YAP Prolog 							 *
+ *									 *
+ *	Yap Prolog was developed at NCCUP - Universidade do Porto	 *
+ *									 *
+ * Copyright L.Damas, V.S.Costa and Universidade do Porto 1985-1997	 *
+ *									 *
+ **************************************************************************
+ *									 *
+ * File:		matrix.c * Last rev:
+ ** mods: * comments:	numerical arrays *
+ *									 *
+ *************************************************************************/
 
-#include "config.h"
 #include "YapInterface.h"
+#include "config.h"
 #include <math.h>
 #if defined(__MINGW32__) || _MSC_VER
 #include <windows.h>
@@ -43,96 +41,82 @@
 /* maximal number of dimensions, 1024 should be enough */
 #define MAX_DIMS 1024
 
-typedef enum {
-  INT_MATRIX,
-  FLOAT_MATRIX
-} mat_data_type;
+typedef enum { INT_MATRIX, FLOAT_MATRIX } mat_data_type;
 
 typedef enum {
-  MAT_TYPE=0,
-  MAT_BASE=1,
-  MAT_NDIMS=2,
-  MAT_SIZE=3,
-  MAT_ALIGN=4,
-  MAT_DIMS=5,
+  MAT_TYPE = 0,
+  MAT_BASE = 1,
+  MAT_NDIMS = 2,
+  MAT_SIZE = 3,
+  MAT_ALIGN = 4,
+  MAT_DIMS = 5,
 } mat_type;
 
 typedef enum {
-  MAT_PLUS=0,
-  MAT_SUB=1,
-  MAT_TIMES=2,
-  MAT_DIV=3,
-  MAT_IDIV=4,
-  MAT_ZDIV=5,
-  MAT_LOG=6,
-  MAT_EXP=7
+  MAT_PLUS = 0,
+  MAT_SUB = 1,
+  MAT_TIMES = 2,
+  MAT_DIV = 3,
+  MAT_IDIV = 4,
+  MAT_ZDIV = 5,
+  MAT_LOG = 6,
+  MAT_EXP = 7
 } op_type;
 
 YAP_Functor FunctorM;
 YAP_Atom AtomC;
 
-static long int *
-matrix_long_data(int *mat, int ndims)
-{
-  return (long int *)(mat+(MAT_DIMS+ndims));
+static long int *matrix_long_data(int *mat, int ndims) {
+  return (long int *)(mat + (MAT_DIMS + ndims));
 }
 
-static double *
-matrix_double_data(int *mat, int ndims)
-{
-  return (double *)(mat+(MAT_DIMS+ndims));
+static double *matrix_double_data(int *mat, int ndims) {
+  return (double *)(mat + (MAT_DIMS + ndims));
 }
 
-static unsigned int
-matrix_get_offset(int *mat, int* indx)
-{
+static unsigned int matrix_get_offset(int *mat, int *indx) {
   unsigned int i, pos = mat[MAT_SIZE], off = 0;
 
   /* find where we are */
   for (i = 0; i < mat[MAT_NDIMS]; i++) {
     int v;
-    pos /= mat[MAT_DIMS+i];
-    v = indx[i]-mat[MAT_BASE];
-    if (v >= mat[MAT_DIMS+i]) {
+    pos /= mat[MAT_DIMS + i];
+    v = indx[i] - mat[MAT_BASE];
+    if (v >= mat[MAT_DIMS + i]) {
       return off;
     }
-    off += pos*v;
+    off += pos * v;
   }
   return off;
 }
 
-static void
-matrix_get_index(int *mat, unsigned int offset, int* indx)
-{
+static void matrix_get_index(int *mat, unsigned int offset, int *indx) {
   unsigned int i, pos = mat[MAT_SIZE];
 
   /* find where we are */
   for (i = 0; i < mat[MAT_NDIMS]; i++) {
-    pos /= mat[MAT_DIMS+i];
+    pos /= mat[MAT_DIMS + i];
     indx[i] = offset / pos;
     offset = offset % pos;
   }
 }
 
-static void
-matrix_next_index(int *dims, int ndims, int* indx)
-{
+static void matrix_next_index(int *dims, int ndims, int *indx) {
   unsigned int i;
 
   /* find where we are */
-  for (i = ndims; i >0; ) {
+  for (i = ndims; i > 0;) {
     i--;
     indx[i]++;
-    if (indx[i]!=dims[i]) return;
+    if (indx[i] != dims[i])
+      return;
     indx[i] = 0;
   }
 }
 
-static YAP_Term
-new_int_matrix(int ndims, int dims[], long int data[])
-{
+static YAP_Term new_int_matrix(int ndims, int dims[], long int data[]) {
   unsigned int sz;
-  unsigned int i, nelems=1;
+  unsigned int i, nelems = 1;
   YAP_Term blob;
   int *mat;
   long int *bdata;
@@ -140,13 +124,15 @@ new_int_matrix(int ndims, int dims[], long int data[])
 
   /* in case we don't have enough room and need to shift the stack, we can't
      really afford to keep a pointer to the global */
-  for (i=0;i< ndims;i++) {
+  for (i = 0; i < ndims; i++) {
     idims[i] = dims[i];
     nelems *= dims[i];
   }
-  sz = ((MAT_DIMS+1)*sizeof(int)+ndims*sizeof(int)+nelems*sizeof(long int))/sizeof(YAP_CELL);
+  sz = ((MAT_DIMS + 1) * sizeof(int) + ndims * sizeof(int) +
+        nelems * sizeof(long int)) /
+       sizeof(YAP_CELL);
   blob = YAP_MkBlobTerm(sz);
-  if (blob == YAP_TermNil())  {
+  if (blob == YAP_TermNil()) {
     return blob;
   }
   mat = (int *)YAP_BlobOfTerm(blob);
@@ -154,20 +140,18 @@ new_int_matrix(int ndims, int dims[], long int data[])
   mat[MAT_BASE] = 0;
   mat[MAT_NDIMS] = ndims;
   mat[MAT_SIZE] = nelems;
-  for (i=0;i< ndims;i++) {
-    mat[MAT_DIMS+i] = idims[i];
+  for (i = 0; i < ndims; i++) {
+    mat[MAT_DIMS + i] = idims[i];
   }
-  bdata = matrix_long_data(mat,ndims);
+  bdata = matrix_long_data(mat, ndims);
   if (data)
-    memmove((void *)bdata,(void *)data,sizeof(double)*nelems);
+    memmove((void *)bdata, (void *)data, sizeof(double) * nelems);
   return blob;
 }
 
-static YAP_Term
-new_float_matrix(int ndims, int dims[], double data[])
-{
-   unsigned int sz;
-  unsigned int i, nelems=1;
+static YAP_Term new_float_matrix(int ndims, int dims[], double data[]) {
+  unsigned int sz;
+  unsigned int i, nelems = 1;
   YAP_Term blob;
   int *mat;
   double *bdata;
@@ -175,11 +159,13 @@ new_float_matrix(int ndims, int dims[], double data[])
 
   /* in case we don't have enough room and need to shift the stack, we can't
      really afford to keep a pointer to the global */
-  for (i=0;i< ndims;i++) {
+  for (i = 0; i < ndims; i++) {
     idims[i] = dims[i];
     nelems *= dims[i];
   }
-  sz = ((MAT_DIMS+1)*sizeof(int)+ndims*sizeof(int)+(nelems+1)*sizeof(double)+(sizeof(YAP_CELL)-1))/sizeof(YAP_CELL);
+  sz = ((MAT_DIMS + 1) * sizeof(int) + ndims * sizeof(int) +
+        (nelems + 1) * sizeof(double) + (sizeof(YAP_CELL) - 1)) /
+       sizeof(YAP_CELL);
   blob = YAP_MkBlobTerm(sz);
   if (blob == YAP_TermNil())
     return blob;
@@ -188,18 +174,16 @@ new_float_matrix(int ndims, int dims[], double data[])
   mat[MAT_BASE] = 0;
   mat[MAT_NDIMS] = ndims;
   mat[MAT_SIZE] = nelems;
-  for (i=0;i< ndims;i++) {
-    mat[MAT_DIMS+i] = idims[i];
+  for (i = 0; i < ndims; i++) {
+    mat[MAT_DIMS + i] = idims[i];
   }
-  bdata = matrix_double_data(mat,ndims);
+  bdata = matrix_double_data(mat, ndims);
   if (data)
-    memmove((void *)bdata,(void *)data,sizeof(double)*nelems);
+    memmove((void *)bdata, (void *)data, sizeof(double) * nelems);
   return blob;
 }
 
-static YAP_Bool
-scan_dims(int ndims, YAP_Term tl, int dims[MAX_DIMS])
-{
+static YAP_Bool scan_dims(int ndims, YAP_Term tl, int dims[MAX_DIMS]) {
   int i;
 
   for (i = 0; i < ndims; i++) {
@@ -229,16 +213,14 @@ scan_dims(int ndims, YAP_Term tl, int dims[MAX_DIMS])
   return TRUE;
 }
 
-static YAP_Bool
-scan_dims_args(int ndims, YAP_Term tl, int dims[MAX_DIMS])
-{
+static YAP_Bool scan_dims_args(int ndims, YAP_Term tl, int dims[MAX_DIMS]) {
   int i;
 
   for (i = 0; i < ndims; i++) {
     YAP_Term th;
     int d;
 
-    th = YAP_ArgOfTerm(2+i, tl);
+    th = YAP_ArgOfTerm(2 + i, tl);
     if (!YAP_IsIntTerm(th)) {
       /* ERROR */
       return FALSE;
@@ -253,9 +235,7 @@ scan_dims_args(int ndims, YAP_Term tl, int dims[MAX_DIMS])
   return TRUE;
 }
 
-static YAP_Bool
-cp_int_matrix(YAP_Term tl,YAP_Term matrix)
-{
+static YAP_Bool cp_int_matrix(YAP_Term tl, YAP_Term matrix) {
   int *mat = (int *)YAP_BlobOfTerm(matrix);
   int i, nelems = mat[MAT_SIZE];
   long int *j = matrix_long_data(mat, mat[MAT_NDIMS]);
@@ -283,9 +263,7 @@ cp_int_matrix(YAP_Term tl,YAP_Term matrix)
   return TRUE;
 }
 
-static YAP_Bool
-cp_float_matrix(YAP_Term tl,YAP_Term matrix)
-{
+static YAP_Bool cp_float_matrix(YAP_Term tl, YAP_Term matrix) {
   int *mat = (int *)YAP_BlobOfTerm(matrix);
   int i, nelems = mat[MAT_SIZE];
   double *j = matrix_double_data(mat, mat[MAT_NDIMS]);
@@ -316,10 +294,7 @@ cp_float_matrix(YAP_Term tl,YAP_Term matrix)
   return TRUE;
 }
 
-
-static YAP_Bool
-set_int_matrix(YAP_Term matrix,long int set)
-{
+static YAP_Bool set_int_matrix(YAP_Term matrix, long int set) {
   int *mat = (int *)YAP_BlobOfTerm(matrix);
   int i, nelems = mat[MAT_SIZE];
   long int *j = matrix_long_data(mat, mat[MAT_NDIMS]);
@@ -330,9 +305,7 @@ set_int_matrix(YAP_Term matrix,long int set)
   return TRUE;
 }
 
-static YAP_Bool
-set_float_matrix(YAP_Term matrix,double set)
-{
+static YAP_Bool set_float_matrix(YAP_Term matrix, double set) {
   int *mat = (int *)YAP_BlobOfTerm(matrix);
   int i, nelems = mat[MAT_SIZE];
   double *j = matrix_double_data(mat, mat[MAT_NDIMS]);
@@ -343,31 +316,27 @@ set_float_matrix(YAP_Term matrix,double set)
   return TRUE;
 }
 
-static YAP_Bool
-new_ints_matrix(void)
-{
+static YAP_Bool new_ints_matrix(void) {
   int ndims = YAP_IntOfTerm(YAP_ARG1);
   YAP_Term tl = YAP_ARG2, out;
   int dims[MAX_DIMS];
-  
+
   if (!scan_dims(ndims, tl, dims))
     return FALSE;
   out = new_int_matrix(ndims, dims, NULL);
   if (out == YAP_TermNil())
     return FALSE;
-  if (!cp_int_matrix(YAP_ARG3,out))
+  if (!cp_int_matrix(YAP_ARG3, out))
     return FALSE;
   return YAP_Unify(YAP_ARG4, out);
 }
 
-static YAP_Bool
-new_ints_matrix_set(void)
-{
+static YAP_Bool new_ints_matrix_set(void) {
   int ndims = YAP_IntOfTerm(YAP_ARG1);
   YAP_Term tl = YAP_ARG2, out, tset = YAP_ARG3;
   int dims[MAX_DIMS];
   long int set;
-  
+
   if (!YAP_IsIntTerm(tset)) {
     return FALSE;
   }
@@ -375,14 +344,12 @@ new_ints_matrix_set(void)
   if (!scan_dims(ndims, tl, dims))
     return FALSE;
   out = new_int_matrix(ndims, dims, NULL);
-  if (!set_int_matrix(out,set))
+  if (!set_int_matrix(out, set))
     return FALSE;
   return YAP_Unify(YAP_ARG4, out);
 }
 
-static YAP_Bool
-new_floats_matrix(void)
-{
+static YAP_Bool new_floats_matrix(void) {
   int ndims = YAP_IntOfTerm(YAP_ARG1);
   YAP_Term tl = YAP_ARG2, out;
   int dims[MAX_DIMS];
@@ -391,19 +358,17 @@ new_floats_matrix(void)
   out = new_float_matrix(ndims, dims, NULL);
   if (out == YAP_TermNil())
     return FALSE;
-  if (!cp_float_matrix(YAP_ARG3,out))
+  if (!cp_float_matrix(YAP_ARG3, out))
     return FALSE;
   return YAP_Unify(YAP_ARG4, out);
 }
 
-static YAP_Bool
-new_floats_matrix_set(void)
-{
+static YAP_Bool new_floats_matrix_set(void) {
   int ndims = YAP_IntOfTerm(YAP_ARG1);
   YAP_Term tl = YAP_ARG2, out, tset = YAP_ARG3;
   int dims[MAX_DIMS];
   double set;
-  
+
   if (!YAP_IsFloatTerm(tset)) {
     return FALSE;
   }
@@ -411,32 +376,29 @@ new_floats_matrix_set(void)
   if (!scan_dims(ndims, tl, dims))
     return FALSE;
   out = new_float_matrix(ndims, dims, NULL);
-  if (!set_float_matrix(out,set))
+  if (!set_float_matrix(out, set))
     return FALSE;
   return YAP_Unify(YAP_ARG4, out);
 }
 
-static YAP_Term
-float_matrix_to_list(int *mat) {
+static YAP_Term float_matrix_to_list(int *mat) {
   double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 
   /* prepare for worst case with double taking two cells */
-  if (YAP_RequiresExtraStack(6*mat[MAT_SIZE])) {
+  if (YAP_RequiresExtraStack(6 * mat[MAT_SIZE])) {
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
     data = matrix_double_data(mat, mat[MAT_NDIMS]);
   }
   return YAP_FloatsToList(data, mat[MAT_SIZE]);
 }
 
-static YAP_Term
-mk_int_list(int nelems, int *data)
-{
+static YAP_Term mk_int_list(int nelems, int *data) {
   YAP_Term tn = YAP_TermNil();
   YAP_Term tf = tn;
   int i = 0;
 
-  for (i = nelems-1; i>= 0; i--) {
-    tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]),tf);
+  for (i = nelems - 1; i >= 0; i--) {
+    tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]), tf);
     if (tf == tn) {
       /* error */
       return tn;
@@ -445,15 +407,13 @@ mk_int_list(int nelems, int *data)
   return tf;
 }
 
-static YAP_Term
-mk_int_list2(int nelems, int base, int *data)
-{
+static YAP_Term mk_int_list2(int nelems, int base, int *data) {
   YAP_Term tn = YAP_TermNil();
   YAP_Term tf = tn;
   int i = 0;
 
-  for (i = nelems-1; i>= 0; i--) {
-    tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]+base),tf);
+  for (i = nelems - 1; i >= 0; i--) {
+    tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i] + base), tf);
     if (tf == tn) {
       /* error */
       return tn;
@@ -462,15 +422,13 @@ mk_int_list2(int nelems, int base, int *data)
   return tf;
 }
 
-static YAP_Term
-mk_rep_int_list(int nelems, int data)
-{
+static YAP_Term mk_rep_int_list(int nelems, int data) {
   YAP_Term tn = YAP_TermNil();
   YAP_Term tf = tn;
   int i = 0;
 
-  for (i = nelems-1; i>= 0; i--) {
-    tf = YAP_MkPairTerm(YAP_MkIntTerm(data),tf);
+  for (i = nelems - 1; i >= 0; i--) {
+    tf = YAP_MkPairTerm(YAP_MkIntTerm(data), tf);
     if (tf == tn) {
       /* error */
       return tn;
@@ -479,15 +437,13 @@ mk_rep_int_list(int nelems, int data)
   return tf;
 }
 
-static YAP_Term
-mk_long_list(int nelems, long int *data)
-{
+static YAP_Term mk_long_list(int nelems, long int *data) {
   YAP_Term tn = YAP_TermNil();
   YAP_Term tf = tn;
   int i = 0;
 
-  for (i = nelems-1; i>= 0; i--) {
-    tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]),tf);
+  for (i = nelems - 1; i >= 0; i--) {
+    tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]), tf);
     if (tf == tn) {
       /* error */
       return tn;
@@ -496,119 +452,97 @@ mk_long_list(int nelems, long int *data)
   return tf;
 }
 
-
-static YAP_Term
-long_matrix_to_list(int *mat) {
+static YAP_Term long_matrix_to_list(int *mat) {
   long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
 
   /* prepare for worst case with longs evrywhere (3cells + 1) */
-  if (YAP_RequiresExtraStack(5*mat[MAT_SIZE])) {
+  if (YAP_RequiresExtraStack(5 * mat[MAT_SIZE])) {
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
     data = matrix_long_data(mat, mat[MAT_NDIMS]);
   }
   return mk_long_list(mat[MAT_SIZE], data);
 }
 
-static YAP_Term
-matrix_access(int *mat, int *indx)
-{
+static YAP_Term matrix_access(int *mat, int *indx) {
   unsigned int off = matrix_get_offset(mat, indx);
-  if (mat[MAT_TYPE]==FLOAT_MATRIX)
-    return YAP_MkFloatTerm((matrix_double_data(mat,mat[MAT_NDIMS]))[off]);
+  if (mat[MAT_TYPE] == FLOAT_MATRIX)
+    return YAP_MkFloatTerm((matrix_double_data(mat, mat[MAT_NDIMS]))[off]);
   else
-    return YAP_MkIntTerm((matrix_long_data(mat,mat[MAT_NDIMS]))[off]);
+    return YAP_MkIntTerm((matrix_long_data(mat, mat[MAT_NDIMS]))[off]);
 }
 
-static void
-matrix_float_set(int *mat, int *indx, double nval)
-{
+static void matrix_float_set(int *mat, int *indx, double nval) {
   unsigned int off = 0;
 
   off = matrix_get_offset(mat, indx);
-  (matrix_double_data(mat,mat[MAT_NDIMS]))[off] = nval;
+  (matrix_double_data(mat, mat[MAT_NDIMS]))[off] = nval;
 }
 
-static void
-matrix_long_set(int *mat, int *indx, long int nval)
-{
+static void matrix_long_set(int *mat, int *indx, long int nval) {
   unsigned int off = matrix_get_offset(mat, indx);
-  (matrix_long_data(mat,mat[MAT_NDIMS]))[off] = nval;
+  (matrix_long_data(mat, mat[MAT_NDIMS]))[off] = nval;
 }
 
-static void
-matrix_float_set_all(int *mat, double nval)
-{
+static void matrix_float_set_all(int *mat, double nval) {
   int i;
-  double *data = matrix_double_data(mat,mat[MAT_NDIMS]);
+  double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 
-  for (i = 0; i< mat[MAT_SIZE]; i++)
+  for (i = 0; i < mat[MAT_SIZE]; i++)
     data[i] = nval;
 }
 
-static void
-matrix_long_set_all(int *mat, long int nval)
-{
+static void matrix_long_set_all(int *mat, long int nval) {
   int i;
-  long int *data = matrix_long_data(mat,mat[MAT_NDIMS]);
+  long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
 
   if (nval == 0) {
-    memset((void *)data,0,sizeof(long int)*mat[MAT_SIZE]);
+    memset((void *)data, 0, sizeof(long int) * mat[MAT_SIZE]);
   } else {
-    for (i = 0; i< mat[MAT_SIZE]; i++)
+    for (i = 0; i < mat[MAT_SIZE]; i++)
       data[i] = nval;
   }
 }
 
-static void
-matrix_float_add(int *mat, int *indx, double nval)
-{
+static void matrix_float_add(int *mat, int *indx, double nval) {
   unsigned int off;
-  double *dat = matrix_double_data(mat,mat[MAT_NDIMS]);
+  double *dat = matrix_double_data(mat, mat[MAT_NDIMS]);
 
   off = matrix_get_offset(mat, indx);
   dat[off] += nval;
 }
 
-static void
-matrix_long_add(int *mat, int *indx, long int nval)
-{
-  long int *dat = matrix_long_data(mat,mat[MAT_NDIMS]);
+static void matrix_long_add(int *mat, int *indx, long int nval) {
+  long int *dat = matrix_long_data(mat, mat[MAT_NDIMS]);
   unsigned int off = matrix_get_offset(mat, indx);
   dat[off] += nval;
 }
 
-static void
-matrix_inc(int *mat, int *indx)
-{
+static void matrix_inc(int *mat, int *indx) {
   unsigned int off = matrix_get_offset(mat, indx);
-  if (mat[MAT_TYPE]==FLOAT_MATRIX)
-    (matrix_double_data(mat,mat[MAT_NDIMS])[off])++;
+  if (mat[MAT_TYPE] == FLOAT_MATRIX)
+    (matrix_double_data(mat, mat[MAT_NDIMS])[off])++;
   else
-    ((matrix_long_data(mat,mat[MAT_NDIMS]))[off])++;
+    ((matrix_long_data(mat, mat[MAT_NDIMS]))[off])++;
 }
 
-static void
-matrix_dec(int *mat, int *indx)
-{
+static void matrix_dec(int *mat, int *indx) {
   unsigned int off = matrix_get_offset(mat, indx);
-  if (mat[MAT_TYPE]==FLOAT_MATRIX)
-    (matrix_double_data(mat,mat[MAT_NDIMS])[off])--;
+  if (mat[MAT_TYPE] == FLOAT_MATRIX)
+    (matrix_double_data(mat, mat[MAT_NDIMS])[off])--;
   else
-    ((matrix_long_data(mat,mat[MAT_NDIMS]))[off])--;
+    ((matrix_long_data(mat, mat[MAT_NDIMS]))[off])--;
 }
 
-static YAP_Term
-matrix_inc2(int *mat, int *indx)
-{
+static YAP_Term matrix_inc2(int *mat, int *indx) {
   unsigned int off = matrix_get_offset(mat, indx);
-  if (mat[MAT_TYPE]==FLOAT_MATRIX) {
-    double *data  = matrix_double_data(mat,mat[MAT_NDIMS]);
+  if (mat[MAT_TYPE] == FLOAT_MATRIX) {
+    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
     double d = data[off];
     d++;
     data[off] = d;
     return YAP_MkFloatTerm(d);
   } else {
-    long int *data  = matrix_long_data(mat,mat[MAT_NDIMS]);
+    long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
     long int d = data[off];
     d++;
     data[off] = d;
@@ -616,18 +550,16 @@ matrix_inc2(int *mat, int *indx)
   }
 }
 
-static YAP_Term
-matrix_dec2(int *mat, int *indx)
-{
+static YAP_Term matrix_dec2(int *mat, int *indx) {
   unsigned int off = matrix_get_offset(mat, indx);
-  if (mat[MAT_TYPE]==FLOAT_MATRIX) {
-    double *data  = matrix_double_data(mat,mat[MAT_NDIMS]);
+  if (mat[MAT_TYPE] == FLOAT_MATRIX) {
+    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
     double d = data[off];
     d--;
     data[off] = d;
     return YAP_MkFloatTerm(d);
   } else {
-    long int *data  = matrix_long_data(mat,mat[MAT_NDIMS]);
+    long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
     long int d = data[off];
     d--;
     data[off] = d;
@@ -635,13 +567,10 @@ matrix_dec2(int *mat, int *indx)
   }
 }
 
-
-static YAP_Bool
-matrix_set(void)
-{
+static YAP_Bool matrix_set(void) {
   int dims[MAX_DIMS], *mat;
   YAP_Term tf;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -674,13 +603,11 @@ matrix_set(void)
   return TRUE;
 }
 
-static YAP_Bool
-matrix_set2(void)
-{
+static YAP_Bool matrix_set2(void) {
   int dims[MAX_DIMS], *mat;
   YAP_Term tf, t = YAP_ARG1;
- 
-  mat = (int *)YAP_BlobOfTerm(YAP_ArgOfTerm(1,t));
+
+  mat = (int *)YAP_BlobOfTerm(YAP_ArgOfTerm(1, t));
   if (!mat) {
     /* Error */
     return FALSE;
@@ -712,12 +639,10 @@ matrix_set2(void)
   return TRUE;
 }
 
-static YAP_Bool
-matrix_set_all(void)
-{
+static YAP_Bool matrix_set_all(void) {
   int *mat;
   YAP_Term tf;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -746,12 +671,10 @@ matrix_set_all(void)
   return TRUE;
 }
 
-static YAP_Bool
-matrix_add(void)
-{
+static YAP_Bool matrix_add(void) {
   int dims[MAX_DIMS], *mat;
   YAP_Term tf;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -784,12 +707,10 @@ matrix_add(void)
   return TRUE;
 }
 
-static YAP_Bool
-do_matrix_access(void)
-{
+static YAP_Bool do_matrix_access(void) {
   int dims[MAX_DIMS], *mat;
   YAP_Term tf;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -803,12 +724,10 @@ do_matrix_access(void)
   return YAP_Unify(tf, YAP_ARG3);
 }
 
-static YAP_Bool
-do_matrix_access2(void)
-{
+static YAP_Bool do_matrix_access2(void) {
   int dims[MAX_DIMS], *mat;
   YAP_Term tf, t = YAP_ARG1;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ArgOfTerm(1, t));
   if (!mat) {
     /* Error */
@@ -822,11 +741,9 @@ do_matrix_access2(void)
   return YAP_Unify(tf, YAP_ARG2);
 }
 
-static YAP_Bool
-do_matrix_inc(void)
-{
+static YAP_Bool do_matrix_inc(void) {
   int dims[MAX_DIMS], *mat;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -840,11 +757,9 @@ do_matrix_inc(void)
   return TRUE;
 }
 
-static YAP_Bool
-do_matrix_dec(void)
-{
+static YAP_Bool do_matrix_dec(void) {
   int dims[MAX_DIMS], *mat;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -858,11 +773,9 @@ do_matrix_dec(void)
   return TRUE;
 }
 
-static YAP_Bool
-do_matrix_inc2(void)
-{
+static YAP_Bool do_matrix_inc2(void) {
   int dims[MAX_DIMS], *mat;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -872,15 +785,12 @@ do_matrix_inc2(void)
     /* Error */
     return FALSE;
   }
-  return
-    YAP_Unify(matrix_inc2(mat, dims), YAP_ARG3);
+  return YAP_Unify(matrix_inc2(mat, dims), YAP_ARG3);
 }
 
-static YAP_Bool
-do_matrix_dec2(void)
-{
+static YAP_Bool do_matrix_dec2(void) {
   int dims[MAX_DIMS], *mat;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -890,16 +800,13 @@ do_matrix_dec2(void)
     /* Error */
     return FALSE;
   }
-  return
-    YAP_Unify(matrix_dec2(mat, dims), YAP_ARG3);
+  return YAP_Unify(matrix_dec2(mat, dims), YAP_ARG3);
 }
 
-static YAP_Bool
-matrix_to_list(void)
-{
+static YAP_Bool matrix_to_list(void) {
   int *mat;
   YAP_Term tf;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -912,11 +819,9 @@ matrix_to_list(void)
   return YAP_Unify(YAP_ARG2, tf);
 }
 
-static YAP_Bool
-matrix_set_base(void)
-{
+static YAP_Bool matrix_set_base(void) {
   int *mat;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -926,43 +831,36 @@ matrix_set_base(void)
   return TRUE;
 }
 
-
-static YAP_Bool
-matrix_dims(void)
-{
+static YAP_Bool matrix_dims(void) {
   int *mat;
   YAP_Term tf;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
     return FALSE;
   }
-  tf = mk_int_list(mat[MAT_NDIMS],mat+MAT_DIMS);
+  tf = mk_int_list(mat[MAT_NDIMS], mat + MAT_DIMS);
   return YAP_Unify(YAP_ARG2, tf);
 }
 
-static YAP_Bool
-matrix_dims3(void)
-{
+static YAP_Bool matrix_dims3(void) {
   int *mat;
   YAP_Term tf, tof;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
     return FALSE;
   }
-  tf = mk_int_list(mat[MAT_NDIMS],mat+MAT_DIMS);
-  tof = mk_rep_int_list(mat[MAT_NDIMS],mat[MAT_BASE]);
+  tf = mk_int_list(mat[MAT_NDIMS], mat + MAT_DIMS);
+  tof = mk_rep_int_list(mat[MAT_NDIMS], mat[MAT_BASE]);
   return YAP_Unify(YAP_ARG2, tf) && YAP_Unify(YAP_ARG3, tof);
 }
 
-static YAP_Bool
-matrix_size(void)
-{
+static YAP_Bool matrix_size(void) {
   int *mat;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -971,11 +869,9 @@ matrix_size(void)
   return YAP_Unify(YAP_ARG2, YAP_MkIntTerm(mat[MAT_SIZE]));
 }
 
-static YAP_Bool
-matrix_ndims(void)
-{
+static YAP_Bool matrix_ndims(void) {
   int *mat;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* Error */
@@ -984,12 +880,10 @@ matrix_ndims(void)
   return YAP_Unify(YAP_ARG2, YAP_MkIntTerm(mat[MAT_NDIMS]));
 }
 
-static YAP_Bool
-matrix_type(void)
-{
+static YAP_Bool matrix_type(void) {
   int *mat;
   YAP_Term tf;
-  
+
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
   if (!mat) {
     /* not an error, it may be called on a term matrix */
@@ -1003,9 +897,7 @@ matrix_type(void)
   return YAP_Unify(YAP_ARG2, tf);
 }
 
-static YAP_Bool
-matrix_arg_to_offset(void)
-{
+static YAP_Bool matrix_arg_to_offset(void) {
   int indx[MAX_DIMS], *mat;
   unsigned int off;
 
@@ -1023,9 +915,7 @@ matrix_arg_to_offset(void)
   return YAP_Unify(YAP_ARG3, YAP_MkIntTerm(off));
 }
 
-static YAP_Bool
-matrix_offset_to_arg(void)
-{
+static YAP_Bool matrix_offset_to_arg(void) {
   int indx[MAX_DIMS], *mat;
   unsigned int off;
   YAP_Term ti, tf;
@@ -1045,65 +935,55 @@ matrix_offset_to_arg(void)
   return YAP_Unify(YAP_ARG3, tf);
 }
 
-static unsigned int
-scan_max_long(int sz, long int *data)
-{
-  int i, off=0;
-  long int max= data[0];
-  for (i=1; i<sz; i++) {
-    if (data[i]>max) {
-      off=i;
+static unsigned int scan_max_long(int sz, long int *data) {
+  int i, off = 0;
+  long int max = data[0];
+  for (i = 1; i < sz; i++) {
+    if (data[i] > max) {
+      off = i;
       max = data[i];
     }
   }
   return off;
 }
 
-static unsigned int
-scan_max_float(int sz, double *data)
-{
-  int i, off=0;
-  double max= data[0];
-  for (i=1; i<sz; i++) {
-    if (data[i]>max) {
+static unsigned int scan_max_float(int sz, double *data) {
+  int i, off = 0;
+  double max = data[0];
+  for (i = 1; i < sz; i++) {
+    if (data[i] > max) {
       max = data[i];
-      off=i;
+      off = i;
     }
   }
   return off;
 }
 
-static unsigned int
-scan_min_long(int sz, long int *data)
-{
-  int i, off=0;
-  long int max= data[0];
-  for (i=1; i<sz; i++) {
-    if (data[i]<max) {
+static unsigned int scan_min_long(int sz, long int *data) {
+  int i, off = 0;
+  long int max = data[0];
+  for (i = 1; i < sz; i++) {
+    if (data[i] < max) {
       max = data[i];
-      off=i;
+      off = i;
     }
   }
   return off;
 }
 
-static unsigned int
-scan_min_float(int sz, double *data)
-{
-  int i, off=0;
-  double max= data[0];
-  for (i=1; i<sz; i++) {
-    if (data[i]<max) {
+static unsigned int scan_min_float(int sz, double *data) {
+  int i, off = 0;
+  double max = data[0];
+  for (i = 1; i < sz; i++) {
+    if (data[i] < max) {
       max = data[i];
-      off=i;
+      off = i;
     }
   }
   return off;
 }
 
-static YAP_Bool
-matrix_max(void)
-{
+static YAP_Bool matrix_max(void) {
   int *mat;
   unsigned int off;
   YAP_Term tf;
@@ -1125,9 +1005,7 @@ matrix_max(void)
   return YAP_Unify(YAP_ARG2, tf);
 }
 
-static YAP_Bool
-matrix_maxarg(void)
-{
+static YAP_Bool matrix_maxarg(void) {
   int indx[MAX_DIMS], *mat;
   unsigned int off;
   YAP_Term tf;
@@ -1149,9 +1027,7 @@ matrix_maxarg(void)
   return YAP_Unify(YAP_ARG2, tf);
 }
 
-static YAP_Bool
-matrix_min(void)
-{
+static YAP_Bool matrix_min(void) {
   int *mat;
   unsigned int off;
   YAP_Term tf;
@@ -1173,9 +1049,7 @@ matrix_min(void)
   return YAP_Unify(YAP_ARG2, tf);
 }
 
-static YAP_Bool
-matrix_log_all(void)
-{
+static YAP_Bool matrix_log_all(void) {
   int *mat;
 
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
@@ -1189,16 +1063,14 @@ matrix_log_all(void)
     double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
     int i;
 
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       data[i] = log(data[i]);
     }
   }
   return TRUE;
 }
 
-static YAP_Bool
-matrix_log_all2(void)
-{
+static YAP_Bool matrix_log_all2(void) {
   int *mat;
 
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
@@ -1216,13 +1088,13 @@ matrix_log_all2(void)
     if (!YAP_IsVarTerm(YAP_ARG2)) {
       out = YAP_ARG2;
     } else {
-      out = new_float_matrix(mat[MAT_NDIMS], mat+MAT_DIMS, NULL);
+      out = new_float_matrix(mat[MAT_NDIMS], mat + MAT_DIMS, NULL);
       if (out == YAP_TermNil())
-	return FALSE;
+        return FALSE;
     }
     nmat = (int *)YAP_BlobOfTerm(out);
     ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       ndata[i] = log((double)data[i]);
     }
     if (YAP_IsVarTerm(YAP_ARG2)) {
@@ -1237,13 +1109,13 @@ matrix_log_all2(void)
     if (!YAP_IsVarTerm(YAP_ARG2)) {
       out = YAP_ARG2;
     } else {
-      out = new_float_matrix(mat[MAT_NDIMS], mat+MAT_DIMS, NULL);
+      out = new_float_matrix(mat[MAT_NDIMS], mat + MAT_DIMS, NULL);
       if (out == YAP_TermNil())
-	return FALSE;
+        return FALSE;
     }
     nmat = (int *)YAP_BlobOfTerm(out);
     ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       ndata[i] = log(data[i]);
     }
     if (YAP_IsVarTerm(YAP_ARG2)) {
@@ -1253,9 +1125,7 @@ matrix_log_all2(void)
   return TRUE;
 }
 
-static YAP_Bool
-matrix_exp_all(void)
-{
+static YAP_Bool matrix_exp_all(void) {
   int *mat;
 
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
@@ -1269,16 +1139,14 @@ matrix_exp_all(void)
     double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
     int i;
 
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       data[i] = exp(data[i]);
     }
   }
   return TRUE;
 }
 
-static YAP_Bool
-matrix_exp2_all(void)
-{
+static YAP_Bool matrix_exp2_all(void) {
   int *mat;
 
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
@@ -1293,20 +1161,19 @@ matrix_exp2_all(void)
     int i;
     double max = data[0];
 
-    for (i=1; i< mat[MAT_SIZE]; i++) {
-      if (data[i] > max) max = data[i];
+    for (i = 1; i < mat[MAT_SIZE]; i++) {
+      if (data[i] > max)
+        max = data[i];
     }
-    
-    for (i=0; i< mat[MAT_SIZE]; i++) {
-      data[i] = exp(data[i]-max);
+
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
+      data[i] = exp(data[i] - max);
     }
   }
   return TRUE;
 }
 
-static YAP_Bool
-matrix_exp_all2(void)
-{
+static YAP_Bool matrix_exp_all2(void) {
   int *mat;
 
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
@@ -1324,13 +1191,13 @@ matrix_exp_all2(void)
     if (!YAP_IsVarTerm(YAP_ARG2)) {
       out = YAP_ARG2;
     } else {
-      out = new_float_matrix(mat[MAT_NDIMS], mat+MAT_DIMS, NULL);
+      out = new_float_matrix(mat[MAT_NDIMS], mat + MAT_DIMS, NULL);
       if (out == YAP_TermNil())
-	return FALSE;
+        return FALSE;
     }
     nmat = (int *)YAP_BlobOfTerm(out);
     ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       ndata[i] = exp((double)data[i]);
     }
     if (YAP_IsVarTerm(YAP_ARG2)) {
@@ -1345,13 +1212,13 @@ matrix_exp_all2(void)
     if (!YAP_IsVarTerm(YAP_ARG2)) {
       out = YAP_ARG2;
     } else {
-      out = new_float_matrix(mat[MAT_NDIMS], mat+MAT_DIMS, NULL);
+      out = new_float_matrix(mat[MAT_NDIMS], mat + MAT_DIMS, NULL);
       if (out == YAP_TermNil())
-	return FALSE;
+        return FALSE;
     }
     nmat = (int *)YAP_BlobOfTerm(out);
     ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       ndata[i] = exp(data[i]);
     }
     if (YAP_IsVarTerm(YAP_ARG2)) {
@@ -1361,9 +1228,7 @@ matrix_exp_all2(void)
   return TRUE;
 }
 
-static YAP_Bool
-matrix_minarg(void)
-{
+static YAP_Bool matrix_minarg(void) {
   int indx[MAX_DIMS], *mat;
   unsigned int off;
   YAP_Term tf;
@@ -1385,9 +1250,7 @@ matrix_minarg(void)
   return YAP_Unify(YAP_ARG2, tf);
 }
 
-static YAP_Bool
-matrix_sum(void)
-{
+static YAP_Bool matrix_sum(void) {
   int *mat;
   YAP_Term tf;
 
@@ -1409,48 +1272,52 @@ matrix_sum(void)
     double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
     int i;
     double sum = 0.0;
-
+    // function KahanSum(input)
+    double c = 0.0; // A running compensation for lost low-order bits.
     for (i = 0; i < mat[MAT_SIZE]; i++) {
-      sum += data[i];
+      double y = data[i] - c; // So far, so good: c is zero.
+      double t =
+          sum +
+          y; // Alas, sum is big, y small, so low-order digits of y are lost.
+      c = (t - sum) - y; // (t - sum) cancels the high-order part of y;
+                         // subtracting y recovers negative (low part of y)
+      sum = t;            // Algebraically, c should always be zero. Beware
+                         // overly-aggressive optimizing compilers!
     }
     tf = YAP_MkFloatTerm(sum);
   }
   return YAP_Unify(YAP_ARG2, tf);
 }
 
-static void
-add_int_lines(int total,int nlines,long int *mat0,long int *matf)
-{
-  int ncols = total/nlines, i;
-  for (i=0;i<ncols;i++) {
+static void add_int_lines(int total, int nlines, long int *mat0,
+                          long int *matf) {
+  int ncols = total / nlines, i;
+  for (i = 0; i < ncols; i++) {
     long int sum = 0;
     int j;
 
-    for (j=i;j<total;j+=ncols) {
+    for (j = i; j < total; j += ncols) {
       sum += mat0[j];
     }
     matf[i] = sum;
   }
 }
 
-static void
-add_double_lines(int total,int nlines,double *mat0,double *matf)
-{
-  int ncols = total/nlines, i;
-  for (i=0;i<ncols;i++) {
+static void add_double_lines(int total, int nlines, double *mat0,
+                             double *matf) {
+  int ncols = total / nlines, i;
+  for (i = 0; i < ncols; i++) {
     double sum = 0;
     int j;
 
-    for (j=i;j<total;j+=ncols) {
+    for (j = i; j < total; j += ncols) {
       sum += mat0[j];
     }
     matf[i] = sum;
   }
 }
 
-static YAP_Bool
-matrix_agg_lines(void)
-{
+static YAP_Bool matrix_agg_lines(void) {
   int *mat;
   YAP_Term tf;
   YAP_Term top = YAP_ARG2;
@@ -1471,15 +1338,15 @@ matrix_agg_lines(void)
     int dims = mat[MAT_NDIMS];
     int *nmat;
 
-    tf = new_int_matrix(dims-1,mat+(MAT_DIMS+1),NULL);
+    tf = new_int_matrix(dims - 1, mat + (MAT_DIMS + 1), NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
     nmat = (int *)YAP_BlobOfTerm(tf);
     data = matrix_long_data(mat, dims);
-    ndata = matrix_long_data(nmat, dims-1);
+    ndata = matrix_long_data(nmat, dims - 1);
     if (op == MAT_PLUS) {
-      add_int_lines(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
+      add_int_lines(mat[MAT_SIZE], mat[MAT_DIMS], data, ndata);
     } else
       return FALSE;
   } else {
@@ -1487,53 +1354,48 @@ matrix_agg_lines(void)
     int dims = mat[MAT_NDIMS];
     int *nmat;
 
-    tf = new_float_matrix(dims-1,mat+(MAT_DIMS+1),NULL);
+    tf = new_float_matrix(dims - 1, mat + (MAT_DIMS + 1), NULL);
     nmat = (int *)YAP_BlobOfTerm(tf);
     if (tf == YAP_TermNil())
       return FALSE;
     data = matrix_double_data(mat, dims);
-    ndata = matrix_double_data(nmat, dims-1);
+    ndata = matrix_double_data(nmat, dims - 1);
     if (op == MAT_PLUS) {
-      add_double_lines(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
+      add_double_lines(mat[MAT_SIZE], mat[MAT_DIMS], data, ndata);
     } else
       return FALSE;
   }
-  return YAP_Unify(YAP_ARG3,tf);
+  return YAP_Unify(YAP_ARG3, tf);
 }
 
-static void
-add_int_cols(int total,int nlines,long int *mat0,long int *matf)
-{
-  int ncols = total/nlines, i, j = 0;
-  for (i=0;i<nlines;i++) {
+static void add_int_cols(int total, int nlines, long int *mat0,
+                         long int *matf) {
+  int ncols = total / nlines, i, j = 0;
+  for (i = 0; i < nlines; i++) {
     long int sum = 0;
-    int max = (i+1)*ncols;
+    int max = (i + 1) * ncols;
 
-    for (;j<max;j++) {
+    for (; j < max; j++) {
       sum += mat0[j];
     }
     matf[i] = sum;
   }
 }
 
-static void
-add_double_cols(int total,int nlines,double *mat0,double *matf)
-{
-  int ncols = total/nlines, i, j = 0;
-  for (i=0;i<nlines;i++) {
+static void add_double_cols(int total, int nlines, double *mat0, double *matf) {
+  int ncols = total / nlines, i, j = 0;
+  for (i = 0; i < nlines; i++) {
     double sum = 0;
-    int max = (i+1)*ncols;
+    int max = (i + 1) * ncols;
 
-    for (;j<max;j++) {
+    for (; j < max; j++) {
       sum += mat0[j];
     }
     matf[i] = sum;
   }
 }
 
-static YAP_Bool
-matrix_agg_cols(void)
-{
+static YAP_Bool matrix_agg_cols(void) {
   int *mat;
   YAP_Term tf;
   YAP_Term top = YAP_ARG2;
@@ -1554,7 +1416,7 @@ matrix_agg_cols(void)
     int dims = mat[MAT_NDIMS];
     int *nmat;
 
-    tf = new_int_matrix(1,mat+MAT_DIMS,NULL);
+    tf = new_int_matrix(1, mat + MAT_DIMS, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
@@ -1562,7 +1424,7 @@ matrix_agg_cols(void)
     data = matrix_long_data(mat, dims);
     ndata = matrix_long_data(nmat, 1);
     if (op == MAT_PLUS) {
-      add_int_cols(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
+      add_int_cols(mat[MAT_SIZE], mat[MAT_DIMS], data, ndata);
     } else
       return FALSE;
   } else {
@@ -1570,59 +1432,53 @@ matrix_agg_cols(void)
     int dims = mat[MAT_NDIMS];
     int *nmat;
 
-    tf = new_float_matrix(1,mat+MAT_DIMS,NULL);
+    tf = new_float_matrix(1, mat + MAT_DIMS, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     nmat = (int *)YAP_BlobOfTerm(tf);
     data = matrix_double_data(mat, dims);
     ndata = matrix_double_data(nmat, 1);
     if (op == MAT_PLUS) {
-      add_double_cols(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
+      add_double_cols(mat[MAT_SIZE], mat[MAT_DIMS], data, ndata);
     } else
       return FALSE;
   }
-  return YAP_Unify(YAP_ARG3,tf);
+  return YAP_Unify(YAP_ARG3, tf);
 }
 
-static void
-div_int_by_lines(int total,int nlines,long int *mat1,long int *mat2,double *ndata)
-{
-  int ncols = total/nlines, i;
-  for (i=0;i<total;i++) {
-    ndata[i] = ((double)mat1[i])/mat2[i%ncols];
+static void div_int_by_lines(int total, int nlines, long int *mat1,
+                             long int *mat2, double *ndata) {
+  int ncols = total / nlines, i;
+  for (i = 0; i < total; i++) {
+    ndata[i] = ((double)mat1[i]) / mat2[i % ncols];
   }
 }
 
-static void
-div_int_by_dlines(int total,int nlines,long int *mat1,double *mat2,double *ndata)
-{
-  int ncols = total/nlines, i;
-  for (i=0;i<total;i++) {
-    ndata[i] = mat1[i]/mat2[i%ncols];
+static void div_int_by_dlines(int total, int nlines, long int *mat1,
+                              double *mat2, double *ndata) {
+  int ncols = total / nlines, i;
+  for (i = 0; i < total; i++) {
+    ndata[i] = mat1[i] / mat2[i % ncols];
   }
 }
 
-static void
-div_float_long_by_lines(int total,int nlines,double *mat1,long int *mat2,double *ndata)
-{
-  int ncols = total/nlines, i;
-  for (i=0;i<total;i++) {
-    ndata[i] = mat1[i]/mat2[i%ncols];
+static void div_float_long_by_lines(int total, int nlines, double *mat1,
+                                    long int *mat2, double *ndata) {
+  int ncols = total / nlines, i;
+  for (i = 0; i < total; i++) {
+    ndata[i] = mat1[i] / mat2[i % ncols];
   }
 }
 
-static void
-div_float_by_lines(int total,int nlines,double *mat1,double *mat2,double *ndata)
-{
-  int ncols = total/nlines, i;
-  for (i=0;i<total;i++) {
-    ndata[i] = mat1[i]/mat2[i%ncols];
+static void div_float_by_lines(int total, int nlines, double *mat1,
+                               double *mat2, double *ndata) {
+  int ncols = total / nlines, i;
+  for (i = 0; i < total; i++) {
+    ndata[i] = mat1[i] / mat2[i % ncols];
   }
 }
 
-static YAP_Bool
-matrix_op_to_lines(void)
-{
+static YAP_Bool matrix_op_to_lines(void) {
   int *mat1, *mat2;
   YAP_Term top = YAP_ARG3;
   op_type op;
@@ -1646,36 +1502,36 @@ matrix_op_to_lines(void)
   if (mat1[MAT_TYPE] == INT_MATRIX) {
     long int *data1;
     int dims = mat1[MAT_NDIMS];
-    int *nmat;   
+    int *nmat;
     data1 = matrix_long_data(mat1, dims);
 
     if (mat2[MAT_TYPE] == INT_MATRIX) {
-      long int *data2 = matrix_long_data(mat2, dims-1);
+      long int *data2 = matrix_long_data(mat2, dims - 1);
       if (op == MAT_DIV) {
-	double *ndata;
+        double *ndata;
 
-	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
-	if (tf == YAP_TermNil())
-	  return FALSE;
-	nmat = YAP_BlobOfTerm(tf);
-	ndata = matrix_double_data(nmat, dims);
-	div_int_by_lines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
+        tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
+        if (tf == YAP_TermNil())
+          return FALSE;
+        nmat = YAP_BlobOfTerm(tf);
+        ndata = matrix_double_data(nmat, dims);
+        div_int_by_lines(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
       } else {
-	return FALSE;
+        return FALSE;
       }
     } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
-      double *data2 = matrix_double_data(mat2, dims-1);
+      double *data2 = matrix_double_data(mat2, dims - 1);
       if (op == MAT_DIV) {
-	double *ndata;
+        double *ndata;
 
-	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
-	if (tf == YAP_TermNil())
-	  return FALSE;
-	nmat = YAP_BlobOfTerm(tf);
-	ndata = matrix_double_data(nmat, dims);
-	div_int_by_dlines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
+        tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
+        if (tf == YAP_TermNil())
+          return FALSE;
+        nmat = YAP_BlobOfTerm(tf);
+        ndata = matrix_double_data(nmat, dims);
+        div_int_by_dlines(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
       } else {
-	return FALSE;
+        return FALSE;
       }
     } else {
       return FALSE;
@@ -1686,229 +1542,209 @@ matrix_op_to_lines(void)
     int *nmat;
 
     data1 = matrix_double_data(mat1, dims);
-    tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL); 
+    tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
     nmat = YAP_BlobOfTerm(tf);
     if (tf == YAP_TermNil())
       return FALSE;
     ndata = matrix_double_data(nmat, dims);
     if (mat2[MAT_TYPE] == INT_MATRIX) {
-      long int *data2 = matrix_long_data(mat2, dims-1);
+      long int *data2 = matrix_long_data(mat2, dims - 1);
       if (op == MAT_DIV) {
-	div_float_long_by_lines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
+        div_float_long_by_lines(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2,
+                                ndata);
       } else {
-	return FALSE;
+        return FALSE;
       }
     } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
-      double *data2 = matrix_double_data(mat2, dims-1);
+      double *data2 = matrix_double_data(mat2, dims - 1);
       if (op == MAT_DIV) {
-	div_float_by_lines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
+        div_float_by_lines(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
       } else {
-	return FALSE;
+        return FALSE;
       }
     } else {
       return FALSE;
     }
   }
-  return YAP_Unify(YAP_ARG4,tf);
+  return YAP_Unify(YAP_ARG4, tf);
 }
 
-
-static void
-matrix_long_add_data(long int *nmat, int siz, long int mat1[], long int mat2[])
-{
+static void matrix_long_add_data(long int *nmat, int siz, long int mat1[],
+                                 long int mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]+mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] + mat2[i];
   }
 }
 
-static void
-matrix_long_double_add_data(double *nmat, int siz, long int mat1[], double mat2[])
-{
+static void matrix_long_double_add_data(double *nmat, int siz, long int mat1[],
+                                        double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]+mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] + mat2[i];
   }
 }
 
-static void
-matrix_double_add_data(double *nmat, int siz, double mat1[], double mat2[])
-{
+static void matrix_double_add_data(double *nmat, int siz, double mat1[],
+                                   double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]+mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] + mat2[i];
   }
 }
 
-static void
-matrix_long_sub_data(long int *nmat, int siz, long int mat1[], long int mat2[])
-{
+static void matrix_long_sub_data(long int *nmat, int siz, long int mat1[],
+                                 long int mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]-mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] - mat2[i];
   }
 }
 
-static void
-matrix_long_double_sub_data(double *nmat, int siz, long int mat1[], double mat2[])
-{
+static void matrix_long_double_sub_data(double *nmat, int siz, long int mat1[],
+                                        double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]-mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] - mat2[i];
   }
 }
 
-static void
-matrix_long_double_rsub_data(double *nmat, int siz, double mat1[], long int mat2[])
-{
+static void matrix_long_double_rsub_data(double *nmat, int siz, double mat1[],
+                                         long int mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat2[i]-mat1[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat2[i] - mat1[i];
   }
 }
 
-static void
-matrix_double_sub_data(double *nmat, int siz, double mat1[], double mat2[])
-{
+static void matrix_double_sub_data(double *nmat, int siz, double mat1[],
+                                   double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]-mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] - mat2[i];
   }
 }
 
-static void
-matrix_long_mult_data(long int *nmat, int siz, long int mat1[], long int mat2[])
-{
+static void matrix_long_mult_data(long int *nmat, int siz, long int mat1[],
+                                  long int mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]*mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] * mat2[i];
   }
 }
 
-static void
-matrix_long_double_mult_data(double *nmat, int siz, long int mat1[], double mat2[])
-{
+static void matrix_long_double_mult_data(double *nmat, int siz, long int mat1[],
+                                         double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]*mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] * mat2[i];
   }
 }
 
-static void
-matrix_double_mult_data(double *nmat, int siz, double mat1[], double mat2[])
-{
+static void matrix_double_mult_data(double *nmat, int siz, double mat1[],
+                                    double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]*mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] * mat2[i];
   }
 }
 
-static void
-matrix_long_div_data(long int *nmat, int siz, long int mat1[], long int mat2[])
-{
+static void matrix_long_div_data(long int *nmat, int siz, long int mat1[],
+                                 long int mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]/mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] / mat2[i];
   }
 }
 
-static void
-matrix_long_double_div_data(double *nmat, int siz, long int mat1[], double mat2[])
-{
+static void matrix_long_double_div_data(double *nmat, int siz, long int mat1[],
+                                        double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]/mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] / mat2[i];
   }
 }
 
-static void
-matrix_long_double_div2_data(double *nmat, int siz, double mat1[], long int mat2[])
-{
+static void matrix_long_double_div2_data(double *nmat, int siz, double mat1[],
+                                         long int mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]/mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] / mat2[i];
   }
 }
 
-static void
-matrix_double_div_data(double *nmat, int siz, double mat1[], double mat2[])
-{
+static void matrix_double_div_data(double *nmat, int siz, double mat1[],
+                                   double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
-    nmat[i] = mat1[i]/mat2[i];
+  for (i = 0; i < siz; i++) {
+    nmat[i] = mat1[i] / mat2[i];
   }
 }
 
-static void
-matrix_long_zdiv_data(long int *nmat, int siz, long int mat1[], long int mat2[])
-{
+static void matrix_long_zdiv_data(long int *nmat, int siz, long int mat1[],
+                                  long int mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
+  for (i = 0; i < siz; i++) {
     if (mat1[i] == 0)
       nmat[i] = 0;
     else
-      nmat[i] = mat1[i]/mat2[i];
+      nmat[i] = mat1[i] / mat2[i];
   }
 }
 
-static void
-matrix_long_double_zdiv_data(double *nmat, int siz, long int mat1[], double mat2[])
-{
+static void matrix_long_double_zdiv_data(double *nmat, int siz, long int mat1[],
+                                         double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
+  for (i = 0; i < siz; i++) {
     if (mat1[i] == 0)
       nmat[i] = 0;
     else
-      nmat[i] = mat1[i]/mat2[i];
+      nmat[i] = mat1[i] / mat2[i];
   }
 }
 
-static void
-matrix_long_double_zdiv2_data(double *nmat, int siz, double mat1[], long int mat2[])
-{
+static void matrix_long_double_zdiv2_data(double *nmat, int siz, double mat1[],
+                                          long int mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
+  for (i = 0; i < siz; i++) {
     if (mat1[i] == 0.0)
       nmat[i] = 0;
     else
-      nmat[i] = mat1[i]/mat2[i];
+      nmat[i] = mat1[i] / mat2[i];
   }
 }
 
-static void
-matrix_double_zdiv_data(double *nmat, int siz, double mat1[], double mat2[])
-{
+static void matrix_double_zdiv_data(double *nmat, int siz, double mat1[],
+                                    double mat2[]) {
   int i;
 
-  for (i=0; i< siz; i++) {
+  for (i = 0; i < siz; i++) {
     if (mat1[i] == 0.0) {
       nmat[i] = 0.0;
     } else {
-      nmat[i] = mat1[i]/mat2[i];
+      nmat[i] = mat1[i] / mat2[i];
     }
   }
 }
 
-static YAP_Bool
-matrix_op(void)
-{
+static YAP_Bool matrix_op(void) {
   int *mat1, *mat2;
   YAP_Term top = YAP_ARG3;
   op_type op;
@@ -1935,7 +1771,7 @@ matrix_op(void)
   if (mat1[MAT_TYPE] == INT_MATRIX) {
     long int *data1;
     int dims = mat1[MAT_NDIMS];
-    int *nmat;   
+    int *nmat;
     data1 = matrix_long_data(mat1, dims);
 
     if (mat2[MAT_TYPE] == INT_MATRIX) {
@@ -1943,72 +1779,72 @@ matrix_op(void)
       long int *ndata;
 
       if (create)
-	tf = new_int_matrix(dims,mat1+MAT_DIMS,NULL);
+        tf = new_int_matrix(dims, mat1 + MAT_DIMS, NULL);
       if (tf == YAP_TermNil()) {
-	return FALSE;
+        return FALSE;
       } else {
-	/* there may have been an overflow */
-	mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
-	data1 = matrix_long_data(mat1, dims);
-	mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
-	data2 = matrix_long_data(mat2, dims);
+        /* there may have been an overflow */
+        mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
+        data1 = matrix_long_data(mat1, dims);
+        mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
+        data2 = matrix_long_data(mat2, dims);
       }
       nmat = YAP_BlobOfTerm(tf);
       ndata = matrix_long_data(nmat, dims);
       switch (op) {
       case MAT_PLUS:
-	matrix_long_add_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_add_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_SUB:
-	matrix_long_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_TIMES:
-	matrix_long_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_DIV:
-	matrix_long_div_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_div_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_ZDIV:
-	matrix_long_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       default:
-	return FALSE;
+        return FALSE;
       }
     } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
       double *data2;
       double *ndata;
 
       if (create)
-	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
+        tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
       if (tf == YAP_TermNil()) {
-	return FALSE;
+        return FALSE;
       } else {
-	/* there may have been an overflow */
-	mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
-	data1 = matrix_long_data(mat1, dims);
-	mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
-	data2  = matrix_double_data(mat2, dims);
+        /* there may have been an overflow */
+        mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
+        data1 = matrix_long_data(mat1, dims);
+        mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
+        data2 = matrix_double_data(mat2, dims);
       }
       nmat = YAP_BlobOfTerm(tf);
       ndata = matrix_double_data(nmat, dims);
       switch (op) {
       case MAT_PLUS:
-	matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_SUB:
-	matrix_long_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_TIMES:
-	matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_DIV:
-	matrix_long_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_ZDIV:
-	matrix_long_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       default:
-	return FALSE;
+        return FALSE;
       }
     } else {
       return FALSE;
@@ -2016,7 +1852,7 @@ matrix_op(void)
   } else {
     double *data1;
     int dims = mat1[MAT_NDIMS];
-    int *nmat;   
+    int *nmat;
     data1 = matrix_double_data(mat1, dims);
 
     if (mat2[MAT_TYPE] == INT_MATRIX) {
@@ -2024,112 +1860,107 @@ matrix_op(void)
       double *ndata;
 
       if (create)
-	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
+        tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
       if (tf == YAP_TermNil()) {
-	return FALSE;
+        return FALSE;
       } else {
-	/* there may have been an overflow */
-	mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
-	data1 = matrix_double_data(mat1, dims);
-	mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
-	data2  = matrix_long_data(mat2, dims);
+        /* there may have been an overflow */
+        mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
+        data1 = matrix_double_data(mat1, dims);
+        mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
+        data2 = matrix_long_data(mat2, dims);
       }
       nmat = YAP_BlobOfTerm(tf);
       ndata = matrix_double_data(nmat, dims);
       switch (op) {
       case MAT_PLUS:
-	matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data2, data1);
-	break;
+        matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data2, data1);
+        break;
       case MAT_SUB:
-	matrix_long_double_rsub_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_double_rsub_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_TIMES:
-	matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data2, data1);
-	break;
+        matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data2, data1);
+        break;
       case MAT_DIV:
-	matrix_long_double_div2_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_double_div2_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_ZDIV:
-	matrix_long_double_zdiv2_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_long_double_zdiv2_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       default:
-	return FALSE;
+        return FALSE;
       }
     } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
       double *data2;
       double *ndata;
 
       if (create)
-	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
+        tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
       if (tf == YAP_TermNil()) {
-	return FALSE;
+        return FALSE;
       } else {
-	/* there may have been an overflow */
-	mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
-	data1 = matrix_double_data(mat1, dims);
-	mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
-	data2  = matrix_double_data(mat2, dims);
+        /* there may have been an overflow */
+        mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
+        data1 = matrix_double_data(mat1, dims);
+        mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
+        data2 = matrix_double_data(mat2, dims);
       }
       nmat = YAP_BlobOfTerm(tf);
       ndata = matrix_double_data(nmat, dims);
       switch (op) {
       case MAT_PLUS:
-	matrix_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_SUB:
-	matrix_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_TIMES:
-	matrix_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_DIV:
-	matrix_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       case MAT_ZDIV:
-	matrix_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
-	break;
+        matrix_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
+        break;
       default:
-	return FALSE;
+        return FALSE;
       }
     } else {
       return FALSE;
     }
   }
-  return YAP_Unify(YAP_ARG4,tf);
+  return YAP_Unify(YAP_ARG4, tf);
 }
 
-static void
-add_int_by_cols(int total,int nlines,long int *mat1,long int *mat2,long int *ndata)
-{
-  int i, ncols = total/nlines;
-  for (i=0;i<total;i++) {
-    ndata[i] = mat1[i] + mat2[i/ncols];
+static void add_int_by_cols(int total, int nlines, long int *mat1,
+                            long int *mat2, long int *ndata) {
+  int i, ncols = total / nlines;
+  for (i = 0; i < total; i++) {
+    ndata[i] = mat1[i] + mat2[i / ncols];
   }
 }
 
-static void
-add_int_by_dcols(int total,int nlines,long int *mat1,double *mat2,double *ndata)
-{
-  int i, ncols = total/nlines;
-  for (i=0;i<total;i++) {
-    ndata[i] = mat1[i] + mat2[i/ncols];
+static void add_int_by_dcols(int total, int nlines, long int *mat1,
+                             double *mat2, double *ndata) {
+  int i, ncols = total / nlines;
+  for (i = 0; i < total; i++) {
+    ndata[i] = mat1[i] + mat2[i / ncols];
   }
 }
 
-static void
-add_double_by_cols(int total,int nlines,double *mat1,double *mat2,double *ndata)
-{
+static void add_double_by_cols(int total, int nlines, double *mat1,
+                               double *mat2, double *ndata) {
   int i;
-  int ncols = total/nlines;
+  int ncols = total / nlines;
 
-  for (i=0;i<total;i++) {
-    ndata[i] = mat1[i] + mat2[i/ncols];
+  for (i = 0; i < total; i++) {
+    ndata[i] = mat1[i] + mat2[i / ncols];
   }
 }
 
-static YAP_Bool
-matrix_op_to_cols(void)
-{
+static YAP_Bool matrix_op_to_cols(void) {
   int *mat1, *mat2;
   YAP_Term top = YAP_ARG3;
   op_type op;
@@ -2152,36 +1983,36 @@ matrix_op_to_cols(void)
   if (mat1[MAT_TYPE] == INT_MATRIX) {
     long int *data1;
     int dims = mat1[MAT_NDIMS];
-    int *nmat;   
+    int *nmat;
     data1 = matrix_long_data(mat1, dims);
 
     if (mat2[MAT_TYPE] == INT_MATRIX) {
       long int *data2 = matrix_long_data(mat2, 1);
       if (op == MAT_PLUS) {
-	long int *ndata;
+        long int *ndata;
 
-	tf = new_int_matrix(dims,mat1+MAT_DIMS,NULL);
-	if (tf == YAP_TermNil())
-	  return FALSE;
-	nmat = YAP_BlobOfTerm(tf);
-	ndata = matrix_long_data(nmat, dims);
-	add_int_by_cols(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
+        tf = new_int_matrix(dims, mat1 + MAT_DIMS, NULL);
+        if (tf == YAP_TermNil())
+          return FALSE;
+        nmat = YAP_BlobOfTerm(tf);
+        ndata = matrix_long_data(nmat, dims);
+        add_int_by_cols(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
       } else {
-	return FALSE;
+        return FALSE;
       }
     } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
       double *data2 = matrix_double_data(mat2, 1);
       if (op == MAT_PLUS) {
-	double *ndata;
+        double *ndata;
 
-	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
-	if (tf == YAP_TermNil())
-	  return FALSE;
-	nmat = YAP_BlobOfTerm(tf);
-	ndata = matrix_double_data(nmat, dims);
-	add_int_by_dcols(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
+        tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
+        if (tf == YAP_TermNil())
+          return FALSE;
+        nmat = YAP_BlobOfTerm(tf);
+        ndata = matrix_double_data(nmat, dims);
+        add_int_by_dcols(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
       } else {
-	return FALSE;
+        return FALSE;
       }
     } else {
       return FALSE;
@@ -2193,7 +2024,7 @@ matrix_op_to_cols(void)
 
     if (mat2[MAT_TYPE] != FLOAT_MATRIX)
       return FALSE;
-    tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
+    tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     nmat = YAP_BlobOfTerm(tf);
@@ -2201,16 +2032,14 @@ matrix_op_to_cols(void)
     data2 = matrix_double_data(mat2, 1);
     ndata = matrix_double_data(nmat, dims);
     if (op == MAT_PLUS) {
-      add_double_by_cols(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
+      add_double_by_cols(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
     } else
       return FALSE;
   }
-  return YAP_Unify(YAP_ARG4,tf);
+  return YAP_Unify(YAP_ARG4, tf);
 }
 
-static YAP_Bool
-matrix_op_to_all(void)
-{
+static YAP_Bool matrix_op_to_all(void) {
   int *mat;
   YAP_Term tf = 0;
   YAP_Term top = YAP_ARG2;
@@ -2228,7 +2057,7 @@ matrix_op_to_all(void)
   }
   if (YAP_IsVarTerm(YAP_ARG4)) {
     create = TRUE;
-  } 
+  }
   /* create a new array with same dimensions */
   if (mat[MAT_TYPE] == INT_MATRIX) {
     long int *data;
@@ -2239,74 +2068,73 @@ matrix_op_to_all(void)
     if (YAP_IsIntTerm(tnum)) {
       long int num;
       long int *ndata;
-      
+
       num = YAP_IntOfTerm(tnum);
       data = matrix_long_data(mat, dims);
       if (create) {
-	tf = new_int_matrix(dims,mat+(MAT_DIMS),NULL);
-	if (tf == YAP_TermNil())
-	  return FALSE;
-	nmat = (int *)YAP_BlobOfTerm(tf);
-	ndata = matrix_long_data(nmat, dims);
+        tf = new_int_matrix(dims, mat + (MAT_DIMS), NULL);
+        if (tf == YAP_TermNil())
+          return FALSE;
+        nmat = (int *)YAP_BlobOfTerm(tf);
+        ndata = matrix_long_data(nmat, dims);
       } else {
-	nmat = mat;
-	ndata = data; 
+        nmat = mat;
+        ndata = data;
       }
       if (op == MAT_PLUS) {
-	int i;
-	
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = data[i] + num;
-	} 
+        int i;
+
+        for (i = 0; i < mat[MAT_SIZE]; i++) {
+          ndata[i] = data[i] + num;
+        }
       } else if (op == MAT_TIMES) {
-	int i;
-	
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = data[i] * num;
-	}
+        int i;
+
+        for (i = 0; i < mat[MAT_SIZE]; i++) {
+          ndata[i] = data[i] * num;
+        }
       } else {
-	return FALSE;
+        return FALSE;
       }
     } else if (YAP_IsFloatTerm(tnum)) {
       double num;
       double *ndata;
-      
-      
+
       num = YAP_FloatOfTerm(tnum);
       if (create) {
-	tf = new_float_matrix(dims,mat+(MAT_DIMS),NULL);
-	if (tf == YAP_TermNil())
-	  return FALSE;
-	nmat = (int *)YAP_BlobOfTerm(tf);
-	ndata = matrix_double_data(nmat, dims);
+        tf = new_float_matrix(dims, mat + (MAT_DIMS), NULL);
+        if (tf == YAP_TermNil())
+          return FALSE;
+        nmat = (int *)YAP_BlobOfTerm(tf);
+        ndata = matrix_double_data(nmat, dims);
       } else {
-	return FALSE;
+        return FALSE;
       }
       data = matrix_long_data(mat, dims);
       if (op == MAT_PLUS) {
-	int i;
-	
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = data[i] + num;
-	}
+        int i;
+
+        for (i = 0; i < mat[MAT_SIZE]; i++) {
+          ndata[i] = data[i] + num;
+        }
       } else if (op == MAT_SUB) {
-	int i;
-	
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = num - data[i];
-	}
+        int i;
+
+        for (i = 0; i < mat[MAT_SIZE]; i++) {
+          ndata[i] = num - data[i];
+        }
       } else if (op == MAT_TIMES) {
-	int i;
-	
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = data[i] * num;
-	}
+        int i;
+
+        for (i = 0; i < mat[MAT_SIZE]; i++) {
+          ndata[i] = data[i] * num;
+        }
       } else if (op == MAT_DIV) {
-	int i;
-	
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = data[i] / num;
-	}
+        int i;
+
+        for (i = 0; i < mat[MAT_SIZE]; i++) {
+          ndata[i] = data[i] / num;
+        }
       }
     } else {
       return FALSE;
@@ -2324,72 +2152,62 @@ matrix_op_to_all(void)
       return FALSE;
     } else {
       if (!create)
-	return FALSE;
+        return FALSE;
       num = (double)YAP_IntOfTerm(tnum);
     }
     data = matrix_double_data(mat, dims);
     if (create) {
-      tf = new_float_matrix(dims,mat+(MAT_DIMS),NULL);
+      tf = new_float_matrix(dims, mat + (MAT_DIMS), NULL);
       if (tf == YAP_TermNil())
-	return FALSE;
+        return FALSE;
       nmat = (int *)YAP_BlobOfTerm(tf);
       ndata = matrix_double_data(nmat, dims);
     } else {
       nmat = mat;
       ndata = data;
     }
-    switch(op) {
-    case MAT_PLUS:
-      {
-	int i;
+    switch (op) {
+    case MAT_PLUS: {
+      int i;
 
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = data[i] + num;
-	}
+      for (i = 0; i < mat[MAT_SIZE]; i++) {
+        ndata[i] = data[i] + num;
       }
-      break;
-     case MAT_SUB:
-      {
-	int i;
+    } break;
+    case MAT_SUB: {
+      int i;
 
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = num-data[i];
-	}
+      for (i = 0; i < mat[MAT_SIZE]; i++) {
+        ndata[i] = num - data[i];
       }
-      break;
-    case MAT_TIMES:
-      {
-	int i;
-	
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = data[i] * num;
-	}
+    } break;
+    case MAT_TIMES: {
+      int i;
+
+      for (i = 0; i < mat[MAT_SIZE]; i++) {
+        ndata[i] = data[i] * num;
       }
-      break;
-    case MAT_DIV:
-      {
-	int i;
-      
-	for (i = 0; i < mat[MAT_SIZE]; i++) {
-	  ndata[i] = data[i] / num;
-	}
-      } 
-      break;
+    } break;
+    case MAT_DIV: {
+      int i;
+
+      for (i = 0; i < mat[MAT_SIZE]; i++) {
+        ndata[i] = data[i] / num;
+      }
+    } break;
     default:
       return FALSE;
     }
   }
   if (create)
-    return YAP_Unify(YAP_ARG4,tf);
-  return YAP_Unify(YAP_ARG4,YAP_ARG1);
+    return YAP_Unify(YAP_ARG4, tf);
+  return YAP_Unify(YAP_ARG4, YAP_ARG1);
 }
 
 /* given a matrix M and a set of dims, build a new reordered matrix to follow
    the new order
-*/ 
-static YAP_Bool
-matrix_transpose(void)
-{
+*/
+static YAP_Bool matrix_transpose(void) {
   int ndims, i, *dims, *dimsn;
   int conv[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
   YAP_Term tconv, tf;
@@ -2401,23 +2219,23 @@ matrix_transpose(void)
   ndims = mat[MAT_NDIMS];
   if (mat[MAT_TYPE] == INT_MATRIX) {
     /* create a new matrix with the same size */
-    tf = new_int_matrix(ndims,mat+MAT_DIMS,NULL);
+    tf = new_int_matrix(ndims, mat + MAT_DIMS, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
   } else {
     /* create a new matrix with the same size */
-    tf = new_float_matrix(ndims,mat+MAT_DIMS,NULL);
+    tf = new_float_matrix(ndims, mat + MAT_DIMS, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
   }
   /* just in case there was an overflow */
   mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-  nmat = (int *)YAP_BlobOfTerm(tf);  
-  dims = mat+MAT_DIMS;
-  dimsn = nmat+MAT_DIMS;
+  nmat = (int *)YAP_BlobOfTerm(tf);
+  dims = mat + MAT_DIMS;
+  dimsn = nmat + MAT_DIMS;
   /* we now have our target matrix, let us grab our conversion matrix */
   tconv = YAP_ARG2;
-  for (i=0; i < ndims; i++) {
+  for (i = 0; i < ndims; i++) {
     YAP_Term th;
     long int j;
 
@@ -2435,31 +2253,31 @@ matrix_transpose(void)
     next is to copy the elements to the new matrix.
   */
   if (mat[MAT_TYPE] == INT_MATRIX) {
-    long int *data = matrix_long_data(mat,ndims);
+    long int *data = matrix_long_data(mat, ndims);
     /* create a new matrix with the same size */
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       long int x = data[i];
       int j;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(mat, i, indx);
-      for (j = 0; j < ndims; j++)  {
-	nindx[j] = indx[conv[j]];
+      for (j = 0; j < ndims; j++) {
+        nindx[j] = indx[conv[j]];
       }
       matrix_long_set(nmat, nindx, x);
     }
   } else {
-    double *data = matrix_double_data(mat,ndims);
+    double *data = matrix_double_data(mat, ndims);
     /* create a new matrix with the same size */
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       double x = data[i];
       long j;
 
       matrix_get_index(mat, i, indx);
-      for (j = 0; j < ndims; j++) 
-	nindx[j] = indx[conv[j]];
+      for (j = 0; j < ndims; j++)
+        nindx[j] = indx[conv[j]];
       matrix_float_set(nmat, nindx, x);
     }
   }
@@ -2468,10 +2286,8 @@ matrix_transpose(void)
 
 /* given a matrix M and a set of dims, fold one of the dimensions of the
    matrix on one of the elements
-*/ 
-static YAP_Bool
-matrix_select(void)
-{
+*/
+static YAP_Bool matrix_select(void) {
   int ndims, i, j, newdims, prdim, leftarg, *dims, indx[MAX_DIMS];
   int nindx[MAX_DIMS];
   YAP_Term tpdim, tdimarg, tf;
@@ -2483,7 +2299,7 @@ matrix_select(void)
   /* we now have our target matrix, let us grab our conversion arguments */
   tpdim = YAP_ARG2;
   ndims = mat[MAT_NDIMS];
-  dims = mat+MAT_DIMS;
+  dims = mat + MAT_DIMS;
   if (!YAP_IsIntTerm(tpdim)) {
     return FALSE;
   }
@@ -2493,42 +2309,42 @@ matrix_select(void)
     return FALSE;
   }
   leftarg = YAP_IntOfTerm(tdimarg);
-  for (i=0, j=0; i< ndims; i++) {
+  for (i = 0, j = 0; i < ndims; i++) {
     if (i != prdim) {
-      nindx[j]= dims[i];
+      nindx[j] = dims[i];
       j++;
     }
   }
-  newdims = ndims-1;
+  newdims = ndims - 1;
   if (mat[MAT_TYPE] == INT_MATRIX) {
     long int *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_int_matrix(newdims,nindx,NULL);
+    tf = new_int_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_long_data(mat,ndims);
-    ndata = matrix_long_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_long_data(mat, ndims);
+    ndata = matrix_long_data(nmat, newdims);
     /* create a new matrix with smaller size */
-    for (i=0; i< nmat[MAT_SIZE]; i++) {
-      int j,k;
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
+      int j, k;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(nmat, i, indx);
-      for (j = 0, k=0; j < newdims; j++,k++)  {
-	if (j == prdim) {
-	  nindx[k] = leftarg;
-	  k++;
-	}
-	nindx[k]= indx[j];
+      for (j = 0, k = 0; j < newdims; j++, k++) {
+        if (j == prdim) {
+          nindx[k] = leftarg;
+          k++;
+        }
+        nindx[k] = indx[j];
       }
       if (k == prdim) {
-	  nindx[k] = leftarg;
+        nindx[k] = leftarg;
       }
       ndata[i] = data[matrix_get_offset(mat, nindx)];
     }
@@ -2536,31 +2352,31 @@ matrix_select(void)
     double *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_float_matrix(newdims,nindx,NULL);
+    tf = new_float_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_double_data(mat,ndims);
-    ndata = matrix_double_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_double_data(mat, ndims);
+    ndata = matrix_double_data(nmat, newdims);
     /* create a new matrix with the same size */
-    for (i=0; i< nmat[MAT_SIZE]; i++) {
-      int j,k;
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
+      int j, k;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(nmat, i, indx);
-      for (j = 0, k=0; j < newdims; j++,k++)  {
-	if (j == prdim) {
-	  nindx[k] = leftarg;
-	  k++;
-	}
-	nindx[k]= indx[j];
+      for (j = 0, k = 0; j < newdims; j++, k++) {
+        if (j == prdim) {
+          nindx[k] = leftarg;
+          k++;
+        }
+        nindx[k] = indx[j];
       }
       if (k == prdim) {
-	  nindx[k] = leftarg;
+        nindx[k] = leftarg;
       }
       ndata[i] = data[matrix_get_offset(mat, nindx)];
     }
@@ -2569,10 +2385,8 @@ matrix_select(void)
 }
 
 /* given a matrix M and a set of N-1 dims, get the first dimension
-*/ 
-static YAP_Bool
-matrix_column(void)
-{
+ */
+static YAP_Bool matrix_column(void) {
   int size, i, ndims, newdims[1];
   int indx[MAX_DIMS];
   int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
@@ -2585,7 +2399,7 @@ matrix_column(void)
   ndims = mat[MAT_NDIMS];
   /* we now have our target matrix, let us grab our conversion arguments */
   tconv = YAP_ARG2;
-  for (i=1; i < ndims; i++) {
+  for (i = 1; i < ndims; i++) {
     YAP_Term th;
 
     if (!YAP_IsPairTerm(tconv))
@@ -2608,29 +2422,29 @@ matrix_column(void)
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_long_data(mat,ndims);
-    ndata = matrix_long_data(nmat,1);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_long_data(mat, ndims);
+    ndata = matrix_long_data(nmat, 1);
     /* create a new matrix with smaller size */
-    for (i=0; i< size; i++) {
-      indx[0]=i;
+    for (i = 0; i < size; i++) {
+      indx[0] = i;
       ndata[i] = data[matrix_get_offset(mat, indx)];
     }
   } else {
     double *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_float_matrix(1,newdims,NULL);
+    tf = new_float_matrix(1, newdims, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_double_data(mat,ndims);
-    ndata = matrix_double_data(nmat,1);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_double_data(mat, ndims);
+    ndata = matrix_double_data(nmat, 1);
     /* create a new matrix with smaller size */
-    for (i=0; i< size; i++) {
-      indx[0]=i;
+    for (i = 0; i < size; i++) {
+      indx[0] = i;
       ndata[i] = data[matrix_get_offset(mat, indx)];
     }
   }
@@ -2638,10 +2452,8 @@ matrix_column(void)
 }
 
 /* given a matrix M and a set of dims, sum out one of the dimensions
-*/ 
-static YAP_Bool
-matrix_sum_out(void)
-{
+ */
+static YAP_Bool matrix_sum_out(void) {
   int ndims, i, j, newdims, prdim;
   int indx[MAX_DIMS], nindx[MAX_DIMS];
   YAP_Term tpdim, tf;
@@ -2657,10 +2469,10 @@ matrix_sum_out(void)
     return FALSE;
   }
   prdim = YAP_IntOfTerm(tpdim);
-  newdims = ndims-1;
-  for (i=0, j=0; i< ndims; i++) {
+  newdims = ndims - 1;
+  for (i = 0, j = 0; i < ndims; i++) {
     if (i != prdim) {
-      nindx[j]= (mat+MAT_DIMS)[i];
+      nindx[j] = (mat + MAT_DIMS)[i];
       j++;
     }
   }
@@ -2668,28 +2480,28 @@ matrix_sum_out(void)
     long int *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_int_matrix(newdims,nindx,NULL);
+    tf = new_int_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_long_data(mat,ndims);
-    ndata = matrix_long_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_long_data(mat, ndims);
+    ndata = matrix_long_data(nmat, newdims);
     /* create a new matrix with smaller size */
-    for (i=0;i<nmat[MAT_SIZE];i++) 
+    for (i = 0; i < nmat[MAT_SIZE]; i++)
       ndata[i] = 0;
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       int j, k;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(mat, i, indx);
-      for (j = 0, k=0; j < ndims; j++)  {
-	if (j != prdim) {
-	  nindx[k++]= indx[j];
-	}
+      for (j = 0, k = 0; j < ndims; j++) {
+        if (j != prdim) {
+          nindx[k++] = indx[j];
+        }
       }
       ndata[matrix_get_offset(nmat, nindx)] += data[i];
     }
@@ -2697,28 +2509,28 @@ matrix_sum_out(void)
     double *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_float_matrix(newdims,nindx,NULL);
+    tf = new_float_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_double_data(mat,ndims);
-    ndata = matrix_double_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_double_data(mat, ndims);
+    ndata = matrix_double_data(nmat, newdims);
     /* create a new matrix with smaller size */
-    for (i=0;i<nmat[MAT_SIZE];i++) 
+    for (i = 0; i < nmat[MAT_SIZE]; i++)
       ndata[i] = 0.0;
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       int j, k;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(mat, i, indx);
-      for (j = 0, k=0; j < ndims; j++)  {
-	if (j != prdim) {
-	  nindx[k++]= indx[j];
-	}
+      for (j = 0, k = 0; j < ndims; j++) {
+        if (j != prdim) {
+          nindx[k++] = indx[j];
+        }
       }
       ndata[matrix_get_offset(nmat, nindx)] += data[i];
     }
@@ -2727,10 +2539,8 @@ matrix_sum_out(void)
 }
 
 /* given a matrix M and a set of dims, sum out one of the dimensions
-*/ 
-static YAP_Bool
-matrix_sum_out_several(void)
-{
+ */
+static YAP_Bool matrix_sum_out_several(void) {
   int ndims, i, *dims, newdims;
   int indx[MAX_DIMS], nindx[MAX_DIMS], conv[MAX_DIMS];
   YAP_Term tf, tconv;
@@ -2740,10 +2550,10 @@ matrix_sum_out_several(void)
     return FALSE;
   }
   ndims = mat[MAT_NDIMS];
-  dims = mat+MAT_DIMS;
+  dims = mat + MAT_DIMS;
   /* we now have our target matrix, let us grab our conversion arguments */
   tconv = YAP_ARG2;
-  for (i=0, newdims=0; i < ndims; i++) {
+  for (i = 0, newdims = 0; i < ndims; i++) {
     YAP_Term th;
 
     if (!YAP_IsPairTerm(tconv))
@@ -2761,69 +2571,69 @@ matrix_sum_out_several(void)
     long int *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_int_matrix(newdims,nindx,NULL);
+    tf = new_int_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_long_data(mat,ndims);
-    ndata = matrix_long_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_long_data(mat, ndims);
+    ndata = matrix_long_data(nmat, newdims);
     /* create a new matrix with smaller size */
-    for (i=0;i<nmat[MAT_SIZE];i++) 
+    for (i = 0; i < nmat[MAT_SIZE]; i++)
       ndata[i] = 0;
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       int j, k;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(mat, i, indx);
-      for (j = 0, k=0; j < ndims; j++)  {
-	if (!conv[j]) {
-	  nindx[k++]= indx[j];
-	}
+      for (j = 0, k = 0; j < ndims; j++) {
+        if (!conv[j]) {
+          nindx[k++] = indx[j];
+        }
       }
-      ndata[matrix_get_offset(nmat, nindx)] = log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
+      ndata[matrix_get_offset(nmat, nindx)] =
+          log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
     }
   } else {
     double *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_float_matrix(newdims,nindx,NULL);
+    tf = new_float_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_double_data(mat,ndims);
-    ndata = matrix_double_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_double_data(mat, ndims);
+    ndata = matrix_double_data(nmat, newdims);
     /* create a new matrix with smaller size */
-    for (i=0;i<nmat[MAT_SIZE];i++) 
+    for (i = 0; i < nmat[MAT_SIZE]; i++)
       ndata[i] = 0.0;
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       int j, k;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(mat, i, indx);
-      for (j = 0, k=0; j < ndims; j++)  {
-	if (!conv[j]) {
-	  nindx[k++]= indx[j];
-	}
+      for (j = 0, k = 0; j < ndims; j++) {
+        if (!conv[j]) {
+          nindx[k++] = indx[j];
+        }
       }
-      ndata[matrix_get_offset(nmat, nindx)] = log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
+      ndata[matrix_get_offset(nmat, nindx)] =
+          log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
     }
   }
   return YAP_Unify(YAP_ARG3, tf);
 }
 
 /* given a matrix M and a set of dims, sum out one of the dimensions
-*/ 
-static YAP_Bool
-matrix_sum_out_logs(void)
-{
+ */
+static YAP_Bool matrix_sum_out_logs(void) {
   int ndims, i, j, *dims, newdims, prdim;
   int nindx[MAX_DIMS];
   YAP_Term tpdim, tf;
@@ -2835,15 +2645,15 @@ matrix_sum_out_logs(void)
   /* we now have our target matrix, let us grab our conversion arguments */
   tpdim = YAP_ARG2;
   ndims = mat[MAT_NDIMS];
-  dims = mat+MAT_DIMS;
+  dims = mat + MAT_DIMS;
   if (!YAP_IsIntTerm(tpdim)) {
     return FALSE;
   }
   prdim = YAP_IntOfTerm(tpdim);
-  newdims = ndims-1;
-  for (i=0, j=0; i< ndims; i++) {
+  newdims = ndims - 1;
+  for (i = 0, j = 0; i < ndims; i++) {
     if (i != prdim) {
-      nindx[j]= dims[i];
+      nindx[j] = dims[i];
       j++;
     }
   }
@@ -2852,28 +2662,28 @@ matrix_sum_out_logs(void)
     int d = 1, j = 0, dd = 1;
 
     /* create a new matrix with the same size */
-    tf = new_int_matrix(newdims,nindx,NULL);
+    tf = new_int_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_long_data(mat,ndims);
-    ndata = matrix_long_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_long_data(mat, ndims);
+    ndata = matrix_long_data(nmat, newdims);
     while (j < prdim) {
-      d = d*dims[j];
+      d = d * dims[j];
       j++;
     }
-    dd = d*dims[prdim];
-    for (i=0;i<nmat[MAT_SIZE];i++) {
-      int j = i % d + (i/dd)*d;
+    dd = d * dims[prdim];
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
+      int j = i % d + (i / dd) * d;
       ndata[j] = exp(data[i]);
     }
-    for (; i< mat[MAT_SIZE]; i++) {
-      int j = i % d + (i/dd)*d;
+    for (; i < mat[MAT_SIZE]; i++) {
+      int j = i % d + (i / dd) * d;
       ndata[j] += exp(data[i]);
     }
-    for (i=0; i< nmat[MAT_SIZE]; i++) {
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
       ndata[i] = log(ndata[i]);
     }
   } else {
@@ -2881,27 +2691,27 @@ matrix_sum_out_logs(void)
     int d = 1, j = 0, dd = 1;
 
     /* create a new matrix with the same size */
-    tf = new_float_matrix(newdims,nindx,NULL);
+    tf = new_float_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_double_data(mat,ndims);
-    ndata = matrix_double_data(nmat,newdims);
-    
-    j = ndims-1;
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_double_data(mat, ndims);
+    ndata = matrix_double_data(nmat, newdims);
+
+    j = ndims - 1;
     while (j > prdim) {
-      d = d*dims[j];
+      d = d * dims[j];
       j--;
     }
-    dd = d*dims[prdim];
-    memset(ndata, 0, sizeof(double)*nmat[MAT_SIZE]);
-    for (i=0; i< mat[MAT_SIZE]; i++) {
-      YAP_Int k = i % d + (i/dd)*d;
+    dd = d * dims[prdim];
+    memset(ndata, 0, sizeof(double) * nmat[MAT_SIZE]);
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
+      YAP_Int k = i % d + (i / dd) * d;
       ndata[k] += exp(data[i]);
     }
-    for (i=0; i< nmat[MAT_SIZE]; i++) {
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
       ndata[i] = log(ndata[i]);
     }
   }
@@ -2909,10 +2719,8 @@ matrix_sum_out_logs(void)
 }
 
 /* given a matrix M and a set of dims, sum out one of the dimensions
-*/ 
-static YAP_Bool
-matrix_sum_out_logs_several(void)
-{
+ */
+static YAP_Bool matrix_sum_out_logs_several(void) {
   int ndims, i, *dims, newdims;
   int indx[MAX_DIMS], nindx[MAX_DIMS], conv[MAX_DIMS];
   YAP_Term tf, tconv;
@@ -2923,10 +2731,10 @@ matrix_sum_out_logs_several(void)
     return FALSE;
   }
   ndims = mat[MAT_NDIMS];
-  dims = mat+MAT_DIMS;
+  dims = mat + MAT_DIMS;
   /* we now have our target matrix, let us grab our conversion arguments */
   tconv = YAP_ARG2;
-  for (i=0, newdims=0; i < ndims; i++) {
+  for (i = 0, newdims = 0; i < ndims; i++) {
     YAP_Term th;
 
     if (!YAP_IsPairTerm(tconv))
@@ -2944,28 +2752,28 @@ matrix_sum_out_logs_several(void)
     long int *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_int_matrix(newdims,nindx,NULL);
+    tf = new_int_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_long_data(mat,ndims);
-    ndata = matrix_long_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_long_data(mat, ndims);
+    ndata = matrix_long_data(nmat, newdims);
     /* create a new matrix with smaller size */
-    for (i=0;i<nmat[MAT_SIZE];i++) 
+    for (i = 0; i < nmat[MAT_SIZE]; i++)
       ndata[i] = 0;
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       int j, k;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(mat, i, indx);
-      for (j = 0, k=0; j < ndims; j++)  {
-	if (!conv[j]) {
-	  nindx[k++]= indx[j];
-	}
+      for (j = 0, k = 0; j < ndims; j++) {
+        if (!conv[j]) {
+          nindx[k++] = indx[j];
+        }
       }
       ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
     }
@@ -2973,32 +2781,32 @@ matrix_sum_out_logs_several(void)
     double *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_float_matrix(newdims,nindx,NULL);
+    tf = new_float_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_double_data(mat,ndims);
-    ndata = matrix_double_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_double_data(mat, ndims);
+    ndata = matrix_double_data(nmat, newdims);
     /* create a new matrix with smaller size */
-    for (i=0;i<nmat[MAT_SIZE];i++) 
+    for (i = 0; i < nmat[MAT_SIZE]; i++)
       ndata[i] = 0.0;
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       int j, k;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(mat, i, indx);
-      for (j = 0, k=0; j < ndims; j++)  {
-	if (!conv[j]) {
-	  nindx[k++]= indx[j];
-	}
+      for (j = 0, k = 0; j < ndims; j++) {
+        if (!conv[j]) {
+          nindx[k++] = indx[j];
+        }
       }
       ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
     }
-    for (i=0; i< nmat[MAT_SIZE]; i++) {
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
       ndata[i] = log(ndata[i]);
     }
   }
@@ -3007,11 +2815,9 @@ matrix_sum_out_logs_several(void)
 
 /* given a matrix M and a set of dims, build a matrix to follow
    the new order
-*/ 
-static YAP_Bool
-matrix_expand(void)
-{
-  int ndims, i, *dims, newdims=0, olddims = 0;
+*/
+static YAP_Bool matrix_expand(void) {
+  int ndims, i, *dims, newdims = 0, olddims = 0;
   int new[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
   YAP_Term tconv, tf;
   int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
@@ -3022,14 +2828,14 @@ matrix_expand(void)
   /* we now have our target matrix, let us grab our conversion matrix */
   tconv = YAP_ARG2;
   ndims = mat[MAT_NDIMS];
-  dims = mat+MAT_DIMS;
-  for (i=0; i < MAX_DIMS; i++) {
+  dims = mat + MAT_DIMS;
+  for (i = 0; i < MAX_DIMS; i++) {
     YAP_Term th;
     long int j;
 
     if (!YAP_IsPairTerm(tconv)) {
       if (tconv != YAP_TermNil())
-	return FALSE;
+        return FALSE;
       break;
     }
     th = YAP_HeadOfTerm(tconv);
@@ -3037,7 +2843,7 @@ matrix_expand(void)
       return FALSE;
     newdims++;
     j = YAP_IntOfTerm(th);
-    if (j==0) {
+    if (j == 0) {
       new[i] = 0;
       nindx[i] = dims[olddims];
       olddims++;
@@ -3053,25 +2859,25 @@ matrix_expand(void)
     long int *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_int_matrix(newdims,nindx,NULL);
+    tf = new_int_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_long_data(mat,ndims);
-    ndata = matrix_long_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_long_data(mat, ndims);
+    ndata = matrix_long_data(nmat, newdims);
     /* create a new matrix with the same size */
-    for (i=0; i< nmat[MAT_SIZE]; i++) {
-      int j,k=0;
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
+      int j, k = 0;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(nmat, i, indx);
-      for (j = 0; j < newdims; j++)  {
-	if (!new[j])
-	  nindx[k++] = indx[j];
+      for (j = 0; j < newdims; j++) {
+        if (!new[j])
+          nindx[k++] = indx[j];
       }
       ndata[i] = data[matrix_get_offset(mat, nindx)];
     }
@@ -3079,29 +2885,29 @@ matrix_expand(void)
     double *data, *ndata;
 
     /* create a new matrix with the same size */
-    tf = new_float_matrix(newdims,nindx,NULL);
+    tf = new_float_matrix(newdims, nindx, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_double_data(mat,ndims);
-    ndata = matrix_double_data(nmat,newdims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_double_data(mat, ndims);
+    ndata = matrix_double_data(nmat, newdims);
     /* create a new matrix with the same size */
-    for (i=0; i < newdims; i++)
+    for (i = 0; i < newdims; i++)
       indx[i] = 0;
-    for (i=0; i< nmat[MAT_SIZE]; i++) {
-      int j,k=0;
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
+      int j, k = 0;
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
-      for (j = 0; j < newdims; j++)  {
-	if (!new[j])
-	  nindx[k++] = indx[j];
+      for (j = 0; j < newdims; j++) {
+        if (!new[j])
+          nindx[k++] = indx[j];
       }
       ndata[i] = data[matrix_get_offset(mat, nindx)];
-      matrix_next_index(nmat+MAT_DIMS, newdims, indx);
+      matrix_next_index(nmat + MAT_DIMS, newdims, indx);
     }
   }
   return YAP_Unify(YAP_ARG3, tf);
@@ -3109,10 +2915,8 @@ matrix_expand(void)
 
 /* given a matrix M and a set of dims, build contract a matrix to follow
    the new order
-*/ 
-static YAP_Bool
-matrix_set_all_that_disagree(void)
-{
+*/
+static YAP_Bool matrix_set_all_that_disagree(void) {
   int ndims, i, *dims;
   int indx[MAX_DIMS];
   YAP_Term tf;
@@ -3125,47 +2929,47 @@ matrix_set_all_that_disagree(void)
     return FALSE;
   }
   ndims = mat[MAT_NDIMS];
-  dims = mat+MAT_DIMS;
+  dims = mat + MAT_DIMS;
   if (mat[MAT_TYPE] == INT_MATRIX) {
     long int *data, *ndata, val;
 
     /* create a new matrix with the same size */
-    tf = new_int_matrix(ndims,dims,NULL);
+    tf = new_int_matrix(ndims, dims, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_long_data(mat,ndims);
-    ndata = matrix_long_data(nmat,ndims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_long_data(mat, ndims);
+    ndata = matrix_long_data(nmat, ndims);
     if (!YAP_IsIntTerm(YAP_ARG4))
       return FALSE;
     val = YAP_IntOfTerm(YAP_ARG4);
     /* create a new matrix with the same size */
-    for (i=0; i< nmat[MAT_SIZE]; i++) {
-      
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
+
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(mat, i, indx);
       if (indx[dim] != pos)
-	ndata[i] = val;
+        ndata[i] = val;
       else
-	ndata[i] = data[i];
+        ndata[i] = data[i];
     }
   } else {
     double *data, *ndata, val;
 
     /* create a new matrix with the same size */
-    tf = new_float_matrix(ndims,dims,NULL);
+    tf = new_float_matrix(ndims, dims, NULL);
     if (tf == YAP_TermNil())
       return FALSE;
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    nmat = (int *)YAP_BlobOfTerm(tf);  
-    data = matrix_double_data(mat,ndims);
-    ndata = matrix_double_data(nmat,ndims);
+    nmat = (int *)YAP_BlobOfTerm(tf);
+    data = matrix_double_data(mat, ndims);
+    ndata = matrix_double_data(nmat, ndims);
     if (YAP_IsFloatTerm(YAP_ARG4))
       val = YAP_FloatOfTerm(YAP_ARG4);
     else if (YAP_IsIntTerm(YAP_ARG4))
@@ -3173,25 +2977,23 @@ matrix_set_all_that_disagree(void)
     else
       return FALSE;
     /* create a new matrix with the same size */
-    for (i=0; i< nmat[MAT_SIZE]; i++) {
-      
+    for (i = 0; i < nmat[MAT_SIZE]; i++) {
+
       /*
-	not very efficient, we could try to take advantage of the fact
-	that we usually only change an index at a time
+        not very efficient, we could try to take advantage of the fact
+        that we usually only change an index at a time
       */
       matrix_get_index(mat, i, indx);
       if (indx[dim] != pos)
-	ndata[i] = val;
+        ndata[i] = val;
       else
-	ndata[i] = data[i];
+        ndata[i] = data[i];
     }
   }
   return YAP_Unify(YAP_ARG5, tf);
 }
 
-static YAP_Bool
-matrix_m(void)
-{
+static YAP_Bool matrix_m(void) {
   int ndims, i, size;
   YAP_Term tm, *tp;
   int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
@@ -3203,10 +3005,10 @@ matrix_m(void)
   size = mat[MAT_SIZE];
   tm = YAP_MkNewApplTerm(FunctorM, 5);
   tp = YAP_ArgsOfTerm(tm);
-  tp[0] = mk_int_list(ndims,mat+MAT_DIMS);
+  tp[0] = mk_int_list(ndims, mat + MAT_DIMS);
   tp[1] = YAP_MkIntTerm(ndims);
   tp[2] = YAP_MkIntTerm(size);
-  tp[3] = mk_rep_int_list(ndims,mat[MAT_BASE]);
+  tp[3] = mk_rep_int_list(ndims, mat[MAT_BASE]);
   tp[4] = YAP_MkNewApplTerm(YAP_MkFunctor(AtomC, size), size);
   tp = YAP_ArgsOfTerm(tp[3]);
   if (mat[MAT_TYPE] == INT_MATRIX) {
@@ -3214,8 +3016,8 @@ matrix_m(void)
 
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    data = matrix_long_data(mat,ndims);
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    data = matrix_long_data(mat, ndims);
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       tp[i] = YAP_MkIntTerm(data[i]);
     }
   } else {
@@ -3223,56 +3025,46 @@ matrix_m(void)
 
     /* in case the matrix moved */
     mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
-    data = matrix_double_data(mat,ndims);
-    for (i=0; i< mat[MAT_SIZE]; i++) {
+    data = matrix_double_data(mat, ndims);
+    for (i = 0; i < mat[MAT_SIZE]; i++) {
       tp[i] = YAP_MkFloatTerm(data[i]);
     }
   }
   return YAP_Unify(YAP_ARG2, tm);
 }
 
-static YAP_Bool
-is_matrix(void)
-{
+static YAP_Bool is_matrix(void) {
   YAP_Term t = YAP_ARG1;
   int *mat = (int *)YAP_BlobOfTerm(t);
 
   if (!mat) {
-    if (!YAP_IsApplTerm(t)) return FALSE;
+    if (!YAP_IsApplTerm(t))
+      return FALSE;
     return YAP_FunctorOfTerm(t) == FunctorM;
   }
   return TRUE;
 }
 
-
-
-static YAP_Bool
-get_float_from_address(void)
-{
+static YAP_Bool get_float_from_address(void) {
   YAP_Float *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
   YAP_Int off = YAP_IntOfTerm(YAP_ARG2);
 
   return YAP_Unify(YAP_ARG3, YAP_MkFloatTerm(fp[off]));
 }
 
-static YAP_Bool
-set_float_from_address(void)
-{
+static YAP_Bool set_float_from_address(void) {
   YAP_Float *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
   YAP_Int off = YAP_IntOfTerm(YAP_ARG2);
-  YAP_Float f = YAP_FloatOfTerm(YAP_ARG3 );
+  YAP_Float f = YAP_FloatOfTerm(YAP_ARG3);
 
   fp[off] = f;
 
   return true;
 }
 
+X_API void init_matrix(void);
 
-X_API void init_matrix( void );
-
-X_API void
-init_matrix(void)
-{
+X_API void init_matrix(void) {
   AtomC = YAP_LookupAtom("c");
   FunctorM = YAP_MkFunctor(YAP_LookupAtom("$matrix"), 5);
 
@@ -3308,16 +3100,18 @@ init_matrix(void)
   YAP_UserCPredicate("matrix_expand", matrix_expand, 3);
   YAP_UserCPredicate("matrix_select", matrix_select, 4);
   YAP_UserCPredicate("matrix_column", matrix_column, 3);
-  YAP_UserCPredicate("matrix_to_logs", matrix_log_all,1);
+  YAP_UserCPredicate("matrix_to_logs", matrix_log_all, 1);
   YAP_UserCPredicate("matrix_to_exps", matrix_exp_all, 1);
   YAP_UserCPredicate("matrix_to_exps2", matrix_exp2_all, 1);
-  YAP_UserCPredicate("matrixn_to_logs", matrix_log_all2,2);
+  YAP_UserCPredicate("matrixn_to_logs", matrix_log_all2, 2);
   YAP_UserCPredicate("matrixn_to_exps", matrix_exp_all2, 2);
   YAP_UserCPredicate("matrix_sum_out", matrix_sum_out, 3);
   YAP_UserCPredicate("matrix_sum_out_several", matrix_sum_out_several, 3);
   YAP_UserCPredicate("matrix_sum_logs_out", matrix_sum_out_logs, 3);
-  YAP_UserCPredicate("matrix_sum_logs_out_several", matrix_sum_out_logs_several, 3);
-  YAP_UserCPredicate("matrix_set_all_that_disagree", matrix_set_all_that_disagree, 5);
+  YAP_UserCPredicate("matrix_sum_logs_out_several", matrix_sum_out_logs_several,
+                     3);
+  YAP_UserCPredicate("matrix_set_all_that_disagree",
+                     matrix_set_all_that_disagree, 5);
   YAP_UserCPredicate("do_matrix_op", matrix_op, 4);
   YAP_UserCPredicate("do_matrix_agg_lines", matrix_agg_lines, 3);
   YAP_UserCPredicate("do_matrix_agg_cols", matrix_agg_cols, 3);
@@ -3326,27 +3120,25 @@ init_matrix(void)
   YAP_UserCPredicate("do_matrix_op_to_cols", matrix_op_to_cols, 4);
   YAP_UserCPredicate("matrix_m", matrix_m, 2);
   YAP_UserCPredicate("matrix", is_matrix, 1);
-  YAP_UserCPredicate("get_float_from_address",get_float_from_address , 3);
-  YAP_UserCPredicate("set_float_from_address",set_float_from_address , 3);
- }
+  YAP_UserCPredicate("get_float_from_address", get_float_from_address, 3);
+  YAP_UserCPredicate("set_float_from_address", set_float_from_address, 3);
+}
 
 #ifdef _WIN32
 
 int WINAPI win_matrixs(HANDLE, DWORD, LPVOID);
 
-int WINAPI win_matrixs(HANDLE hinst, DWORD reason, LPVOID reserved)
-{
-  switch (reason) 
-    {
-    case DLL_PROCESS_ATTACH:
-      break;
-    case DLL_PROCESS_DETACH:
-      break;
-    case DLL_THREAD_ATTACH:
-      break;
-    case DLL_THREAD_DETACH:
-      break;
-    }
+int WINAPI win_matrixs(HANDLE hinst, DWORD reason, LPVOID reserved) {
+  switch (reason) {
+  case DLL_PROCESS_ATTACH:
+    break;
+  case DLL_PROCESS_DETACH:
+    break;
+  case DLL_THREAD_ATTACH:
+    break;
+  case DLL_THREAD_DETACH:
+    break;
+  }
   return 1;
 }
 #endif
diff --git a/packages/python/yap_kernel/yap_ipython/yapi.py b/packages/python/yap_kernel/yap_ipython/yapi.py
index 288564053..d8e4eae1e 100644
--- a/packages/python/yap_kernel/yap_ipython/yapi.py
+++ b/packages/python/yap_kernel/yap_ipython/yapi.py
@@ -116,7 +116,7 @@ class YAPInputSplitter(InputSplitter):
             line = text.rstrip()
         self.errors = []
         engine.mgoal(errors(self, line),"user",True)
-        return self.errors != []
+        return self.errors == []
 
 
     def reset(self):
@@ -496,12 +496,13 @@ class YAPCompleter(Completer):
             completions are coming from different sources this function does not
             ensure that each completion object will only be present once.
         """
+        if text[0] == '%' and text.find('\n') > offset:
+            magic_res = self.magic_matches(text)
+            return text,  magic_res
         self.matches = []
         prolog_res = self.shell.yapeng.mgoal(completions(text, self), "user",True)
-        if self.matches:
-            return text, self.matches
-        magic_res = self.magic_matches(text)
-        return text,  magic_res
+        return text, self.matches
+        
 
 
 
@@ -643,6 +644,7 @@ class YAPRun:
         # variable names should match strings
         # ask = True
         # launch the query
+        cell = raw_cell  # cell has to exist so it can be stored/logged
 
 
         info = interactiveshell.ExecutionInfo(
@@ -679,7 +681,6 @@ class YAPRun:
         #     cell = self.shell.input_transformer_manager.transform_cell(raw_cell)
         # except SyntaxError:
         #     preprocessing_exc_tuple = self.shell.syntax_error() # sys.exc_info()
-        cell = raw_cell  # cell has to exist so it can be stored/logged
         for i in self.errors:
             try:
                 (_,lin,pos,text) = i
@@ -707,31 +708,23 @@ class YAPRun:
         # compiler = self.shell.compile if shell_futures else CachingCompiler()
         self.cell_name = str( self.shell.execution_count)
         if cell[0] == '%':
+            txt0 = cell.split(maxsplit = 1, sep = '\n')
             if cell[1] == '%':
                 linec = False
-                mcell = cell.lstrip('%%')
+                magic = txt[0].lstrip('%%').strip()
+                body = txt[1]
+                self.result = True, self.shell.run_cell_magic(magic, line, body)
             else:
                 linec = True
-                mcell = cell.lstrip('%')
-            txt0 = mcell.split(maxsplit = 2, sep = '\n')
+                magic = cell.lstrip('%').strip()
             txt = txt0[0].split(maxsplit = 2)
             magic = txt[0]
-            if len(txt) == 2:
-                line = txt[1]
-            else:
-                line = ""
-            if linec:
-                self.shell.run_line_magic(magic, line)
-            else:
-                if len(txt0) == 1:
-                    cell = ""
-                else:
-                    body = txt0[1]+'\n'+txt0[2]
-                self.result = True, self.shell.run_cell_magic(magic, line, body)
-                return self.result
+            self.shell.run_line_magic(magic, line)
         # Give the displayhook a reference to our ExecutionResult so it
         # can fill in the output value.
         self.shell.displayhook.exec_result = self.result
+        if syntaxErrors(self, text):
+            self.result.result = False
         has_raised = False
         try:
             self.yapeng.mgoal(streams(True),"user", True)
@@ -743,7 +736,7 @@ class YAPRun:
                 #     trace=1,
                 #     count=0)
                 #
-
+                
                 # def f(self, cell, state):
                 #     state = self.jupyter_query( cell )