diff -p -up iml-1.0.2/src/nullspace.c.orig iml-1.0.2/src/nullspace.c
--- iml-1.0.2/src/nullspace.c.orig	2009-05-08 20:11:29.000000000 -0300
+++ iml-1.0.2/src/nullspace.c	2009-05-08 20:11:53.000000000 -0300
@@ -187,3 +187,150 @@ nullspaceLong(const long n, const long m
     return s;
 
 }
+
+/*
+ * Calling Sequence:
+ *   nullspaceLong(n, m, A, mp_N_pass)
+ *
+ * Summary: Compute the right nullspace of A. In this function A is a
+ * 1-dimensional mpz_t array.
+ *
+ * Input:  n: long, row dimension of A
+ *         m: long, column dimension of A
+ *         A: 1-dim mpz_t array length n*m, representing n x m matrix
+ *            in row major order
+ *
+ * Output:
+ *   - *mp_N_pass: points to a 1-dim mpz_t array of length m*s, where s is the 
+ *                dimension of the right nullspace of A
+ *   - the dimension s of the nullspace is returned
+ *
+ * Notes:
+ *   - The matrix A is represented by one-dimension array in row major order.
+ *   - Space for what mp_N_points to is allocated by this procedure: if the
+ *     nullspace is empty, mp_N_pass is set to NULL.
+ */
+
+long
+nullspaceMP(const long n, const long m, const mpz_t *A, mpz_t * *mp_N_pass)
+{
+    long i, j, k, r, s, *P, *rp, *Pt, *rpt, flag, temp;
+    double *DA;
+    FiniteField p, d = 1;
+    mpz_t *mp_B, *mp_N, mp_D, mp_t1, mp_t2, *C, mp_r;
+
+    mpz_init(mp_r);
+
+    P = XCALLOC(long, n + 1);
+    rp = XCALLOC(long, n + 1);
+    while (1) {
+	p = RandPrime(15, 19);
+	DA = XCALLOC(double, n * m);
+	for (i = 0; i < n * m; i++) {
+  	    mpz_mod_ui (mp_r, A[i], p);
+	    DA[i] = mpz_get_d(mp_r);
+	}
+	for (i = 0; i < n + 1; i++) {
+	    P[i] = i;
+	    rp[i] = 0;
+	}
+	d = 1;
+	RowEchelonTransform(p, DA, n, m, 1, 1, 0, 0, P, rp, &d);
+	XFREE(DA);
+	r = rp[0];
+	s = m - r;
+	if (s == 0) {
+	    *mp_N_pass = NULL;
+	} else if (r == 0) {
+	    flag = 1;
+	    for (i = 0; i < n * m; i++)
+		if ( mpz_cmp_si(A[i],0) )
+		    flag = 0;
+	    if (!flag)
+		continue;
+	    mp_N = XCALLOC(mpz_t, m * m);
+	    for (i = 0; i < m; i++) {
+		for (j = 0; j < m; j++)
+		    mpz_init_set_ui(mp_N[i * m + j], 0);
+		mpz_set_ui(mp_N[i * m + i], 1);
+	    }
+	    *mp_N_pass = mp_N;
+	} else {		/* r>0 and s>0 */
+
+	    Pt = revseq(r, n, P);
+	    rpt = revseq(r, m, rp);
+
+	    C = XCALLOC(mpz_t, r * r);
+	    for (i = 0; i < r; i++)
+		for (j = 0; j < r; j++)
+		    mpz_init_set(C[i * r + j], A[Pt[i] * m + rpt[j]]);
+
+	    mp_B = XCALLOC(mpz_t, r * s);
+	    for (i = 0; i < r; i++)
+		for (j = 0; j < s; j++)
+		    mpz_init_set(mp_B[i * s + j], A[Pt[i] * m + rpt[r + j]]);
+
+	    mpz_init(mp_D);
+	    mp_N = XCALLOC(mpz_t, m * s);
+	    for (i = 0; i < m * s; i++)
+		mpz_init(mp_N[i]);
+
+
+	    nonsingSolvLlhsMM(RightSolu, r, s, C, mp_B, mp_N, mp_D);
+
+	    mpz_neg(mp_D, mp_D);
+	    for (i = 0; i < s; i++)
+		mpz_set(mp_N[(r + i) * s + i], mp_D);
+
+	    for (i = 0; i < r*r; i++)
+	        mpz_clear(C[i]);
+	    XFREE(C);
+
+	    for (i = 0; i < r * s; i++)
+		mpz_clear(mp_B[i]);
+	    XFREE(mp_B);
+	    mpz_clear(mp_D);
+
+	    for (i = r; i >= 1; i--)
+		for (j = 0; j < s; j++)
+		    mpz_swap(mp_N[(i - 1) * s + j],
+			     mp_N[(rp[i] - 1) * s + j]);
+
+	    *mp_N_pass = mp_N;
+
+	    flag = 1;
+	    mpz_init(mp_t1);
+	    mpz_init(mp_t2);
+	    for (i = r; i < n && flag; i++) {
+		for (j = 0; j < s && flag; j++) {
+		    mpz_set_ui(mp_t2, 0);
+		    for (k = 0; k < m; k++) {
+			mpz_mul(mp_t1, mp_N[k * s + j],  A[Pt[i] * m + k]);
+			mpz_add(mp_t2, mp_t2, mp_t1);
+		    }
+		    if (mpz_sgn(mp_t2) != 0)
+			flag = 0;
+		}
+	    }
+	    mpz_clear(mp_t1);
+	    mpz_clear(mp_t2);
+
+	    XFREE(Pt);
+	    XFREE(rpt);
+
+	    if (!flag) {
+		for (i = 0; i < m * s; i++)
+		    mpz_clear(mp_N[i]);
+		XFREE(mp_N);
+		continue;
+	    }
+	}
+	break;
+    }
+    XFREE(P);
+    XFREE(rp);
+
+    mpz_clear(mp_r);
+    return s;
+
+}
diff -p -up iml-1.0.2/src/nullspace.h.orig iml-1.0.2/src/nullspace.h
--- iml-1.0.2/src/nullspace.h.orig	2009-05-08 20:12:44.000000000 -0300
+++ iml-1.0.2/src/nullspace.h	2009-05-08 20:13:10.000000000 -0300
@@ -58,5 +58,9 @@ long nullspaceLong(const long n, 
 		   const long *A, 
 		   mpz_t * *mp_N_pass);
 
+long nullspaceMP(const long n,
+		 const long m,
+		 const mpz_t *A,
+		 mpz_t * *mp_N_pass);
 
 #endif /* !IML_NULLSPACE_H */