diff --git a/Makefile.am b/Makefile.am
index 93d55fff..e5adddf0 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -3,6 +3,13 @@ check_PROGRAMS = mktest
 EXTRA_PROGRAMS = mktables
 
 snapraid_SOURCES = \
+	raid/raid.c \
+	raid/module.c \
+	raid/tables.c \
+	raid/int.c \
+	raid/x86.c \
+	raid/intz.c \
+	raid/x86z.c \
 	snapraid.c \
 	tommy.c \
 	util.c \
@@ -21,7 +28,6 @@ snapraid_SOURCES = \
 	parity.c \
 	handle.c \
 	fnmatch.c \
-	raid.c \
 	selftest.c \
 	speed.c \
 	import.c \
@@ -29,9 +35,13 @@ snapraid_SOURCES = \
 	unix.c
 
 mktables_SOURCES = \
-	mktables.c
+	raid/mktables.c
 
 noinst_HEADERS = \
+	raid/raid.h \
+	raid/internal.h \
+	raid/cpu.h \
+	raid/gf.h \
 	portable.h \
 	snapraid.h \
 	util.h \
@@ -58,10 +68,7 @@ noinst_HEADERS = \
 	tommyhashdyn.c \
 	tommylist.c \
 	fnmatch.h \
-	cpu.h \
-	raid.h \
 	combo.h \
-	tables.h \
 	import.h \
 	mingw.h \
 	unix.h
diff --git a/check.c b/check.c
index 1047314f..ad4c62a2 100644
--- a/check.c
+++ b/check.c
@@ -23,7 +23,7 @@
 #include "state.h"
 #include "parity.h"
 #include "handle.h"
-#include "raid.h"
+#include "raid/raid.h"
 #include "combo.h"
 
 /****************************************************************************/
diff --git a/cpu.h b/raid/cpu.h
similarity index 99%
rename from cpu.h
rename to raid/cpu.h
index 4d8c2079..62504820 100644
--- a/cpu.h
+++ b/raid/cpu.h
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2011,2013 Andrea Mazzoleni
+ * Copyright (C) 2013 Andrea Mazzoleni
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
diff --git a/raid/gf.h b/raid/gf.h
new file mode 100644
index 00000000..2eaedfe5
--- /dev/null
+++ b/raid/gf.h
@@ -0,0 +1,144 @@
+/*
+ * Copyright (C) 2013 Andrea Mazzoleni
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef __RAID_GF_H
+#define __RAID_GF_H
+
+/*
+ * Galois field operations.
+ *
+ * Basic range checks are implemented using BUG_ON().
+ */
+
+/*
+ * Generic assert() with inverted meaning.
+ */
+#define BUG_ON(a) assert(!(a))
+
+
+/*
+ * GF a*b.
+ */
+static __always_inline uint8_t mul(uint8_t a, uint8_t b)
+{
+	return gfmul[a][b];
+}
+
+/*
+ * GF 1/a.
+ * Not defined for a == 0.
+ */
+static __always_inline uint8_t inv(uint8_t v)
+{
+	BUG_ON(v == 0); /* division by zero */
+
+	return gfinv[v];
+}
+
+/*
+ * GF 2^a.
+ */
+static __always_inline uint8_t pow2(int v)
+{
+	BUG_ON(v < 0 || v > 254); /* invalid exponent */
+
+	return gfexp[v];
+}
+
+/*
+ * Gets the multiplication table for a specified value.
+ */
+static __always_inline const uint8_t* table(uint8_t v)
+{
+	return gfmul[v];
+}
+
+extern const uint8_t (*raid_gen)[256];
+
+/*
+ * Gets the generator matrix coefficient for parity 'p' and disk 'd'.
+ */
+static __always_inline uint8_t A(int p, int d)
+{
+	return raid_gen[p][d];
+}
+
+/*
+ * Dereference as uint8_t
+ */
+#define v_8(p) (*(uint8_t*)&(p))
+
+/*
+ * Dereference as uint32_t
+ */
+#define v_32(p) (*(uint32_t*)&(p))
+
+/*
+ * Dereference as uint64_t
+ */
+#define v_64(p) (*(uint64_t*)&(p))
+
+/*
+ * Multiply each byte of a uint32 by 2 in the GF(2^8).
+ */
+static __always_inline uint32_t x2_32(uint32_t v)
+{
+	uint32_t mask = v & 0x80808080U;
+	mask = (mask << 1) - (mask >> 7);
+	v = (v << 1) & 0xfefefefeU;
+	v ^= mask & 0x1d1d1d1dU;
+	return v;
+}
+
+/*
+ * Multiply each byte of a uint64 by 2 in the GF(2^8).
+ */
+static __always_inline uint64_t x2_64(uint64_t v)
+{
+	uint64_t mask = v & 0x8080808080808080ULL;
+	mask = (mask << 1) - (mask >> 7);
+	v = (v << 1) & 0xfefefefefefefefeULL;
+	v ^= mask & 0x1d1d1d1d1d1d1d1dULL;
+	return v;
+}
+
+/*
+ * Divide each byte of a uint32 by 2 in the GF(2^8).
+ */
+static __always_inline uint32_t d2_32(uint32_t v)
+{
+	uint32_t mask = v & 0x01010101U;
+	mask = (mask << 8) - mask;
+	v = (v >> 1) & 0x7f7f7f7fU;
+	v ^= mask & 0x8e8e8e8eU;
+	return v;
+}
+
+/*
+ * Divide each byte of a uint64 by 2 in the GF(2^8).
+ */
+static __always_inline uint64_t d2_64(uint64_t v)
+{
+	uint64_t mask = v & 0x0101010101010101ULL;
+	mask = (mask << 8) - mask;
+	v = (v >> 1) & 0x7f7f7f7f7f7f7f7fULL;
+	v ^= mask & 0x8e8e8e8e8e8e8e8eULL;
+	return v;
+}
+
+#endif
+
diff --git a/raid/int.c b/raid/int.c
new file mode 100644
index 00000000..21417021
--- /dev/null
+++ b/raid/int.c
@@ -0,0 +1,570 @@
+/*
+ * Copyright (C) 2013 Andrea Mazzoleni
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "internal.h"
+#include "gf.h"
+
+/*
+ * PAR1 (RAID5 with xor) 32bit C implementation
+ */
+void raid_par1_int32(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	int d, l;
+	size_t i;
+
+	uint32_t p0;
+	uint32_t p1;
+
+	l = nd - 1;
+	p = v[nd];
+
+	for(i=0;i<size;i+=8) {
+		p0 = v_32(v[l][i]);
+		p1 = v_32(v[l][i+4]);
+		/* accessing disks in backward order because the buffers */
+		/* are also in backward order */
+		for(d=l-1;d>=0;--d) {
+			p0 ^= v_32(v[d][i]);
+			p1 ^= v_32(v[d][i+4]);
+		}
+		v_32(p[i]) = p0;
+		v_32(p[i+4]) = p1;
+	}
+}
+
+/*
+ * PAR1 (RAID5 with xor) 64bit C implementation
+ */
+void raid_par1_int64(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	int d, l;
+	size_t i;
+
+	uint64_t p0;
+	uint64_t p1;
+
+	l = nd - 1;
+	p = v[nd];
+
+	for(i=0;i<size;i+=16) {
+		p0 = v_64(v[l][i]);
+		p1 = v_64(v[l][i+8]);
+		/* accessing disks in backward order because the buffers */
+		/* are also in backward order */
+		for(d=l-1;d>=0;--d) {
+			p0 ^= v_64(v[d][i]);
+			p1 ^= v_64(v[d][i+8]);
+		}
+		v_64(p[i]) = p0;
+		v_64(p[i+8]) = p1;
+	}
+}
+
+/*
+ * PAR2 (RAID6 with powers of 2) 32bit C implementation
+ */
+void raid_par2_int32(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	int d, l;
+	size_t i;
+
+	uint32_t d0, q0, p0;
+	uint32_t d1, q1, p1;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+
+	for(i=0;i<size;i+=8) {
+		q0 = p0 = v_32(v[l][i]);
+		q1 = p1 = v_32(v[l][i+4]);
+		for(d=l-1;d>=0;--d) {
+			d0 = v_32(v[d][i]);
+			d1 = v_32(v[d][i+4]);
+
+			p0 ^= d0;
+			p1 ^= d1;
+
+			q0 = x2_32(q0);
+			q1 = x2_32(q1);
+
+			q0 ^= d0;
+			q1 ^= d1;
+		}
+		v_32(p[i]) = p0;
+		v_32(p[i+4]) = p1;
+		v_32(q[i]) = q0;
+		v_32(q[i+4]) = q1;
+	}
+}
+
+/*
+ * PAR2 (RAID6 with powers of 2) 64bit C implementation
+ */
+void raid_par2_int64(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	int d, l;
+	size_t i;
+
+	uint64_t d0, q0, p0;
+	uint64_t d1, q1, p1;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+
+	for(i=0;i<size;i+=16) {
+		q0 = p0 = v_64(v[l][i]);
+		q1 = p1 = v_64(v[l][i+8]);
+		for(d=l-1;d>=0;--d) {
+			d0 = v_64(v[d][i]);
+			d1 = v_64(v[d][i+8]);
+
+			p0 ^= d0;
+			p1 ^= d1;
+
+			q0 = x2_64(q0);
+			q1 = x2_64(q1);
+
+			q0 ^= d0;
+			q1 ^= d1;
+		}
+		v_64(p[i]) = p0;
+		v_64(p[i+8]) = p1;
+		v_64(q[i]) = q0;
+		v_64(q[i+8]) = q1;
+	}
+}
+
+/*
+ * PAR3 (triple parity with Cauchy matrix) 8bit C implementation
+ * 
+ * Note that instead of a generic multiplicationt table, likely resulting
+ * in multiple cache misses, a precomputed table could be used. 
+ * But this is only a kind of reference function, and we are not really
+ * interested in speed.
+ */
+void raid_par3_int8(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	uint8_t* r;
+	int d, l;
+	size_t i;
+
+	uint8_t d0, r0, q0, p0;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+	r = v[nd+2];
+
+	for(i=0;i<size;i+=1) {
+		p0 = q0 = r0 = 0;
+		for(d=l;d>0;--d) {
+			d0 = v_8(v[d][i]);
+
+			p0 ^= d0;
+			q0 ^= gfmul[d0][gfcauchy[1][d]];
+			r0 ^= gfmul[d0][gfcauchy[2][d]];
+		}
+
+		/* first disk with all coefficients at 1 */
+		d0 = v_8(v[0][i]);
+
+		p0 ^= d0;
+		q0 ^= d0;
+		r0 ^= d0;
+
+		v_8(p[i]) = p0;
+		v_8(q[i]) = q0;
+		v_8(r[i]) = r0;
+	}
+}
+
+/*
+ * PAR4 (quad parity with Cauchy matrix) 8bit C implementation
+ * 
+ * Note that instead of a generic multiplicationt table, likely resulting
+ * in multiple cache misses, a precomputed table could be used. 
+ * But this is only a kind of reference function, and we are not really
+ * interested in speed.
+ */
+void raid_par4_int8(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	uint8_t* r;
+	uint8_t* s;
+	int d, l;
+	size_t i;
+
+	uint8_t d0, s0, r0, q0, p0;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+	r = v[nd+2];
+	s = v[nd+3];
+
+	for(i=0;i<size;i+=1) {
+		p0 = q0 = r0 = s0 = 0;
+		for(d=l;d>0;--d) {
+			d0 = v_8(v[d][i]);
+
+			p0 ^= d0;
+			q0 ^= gfmul[d0][gfcauchy[1][d]];
+			r0 ^= gfmul[d0][gfcauchy[2][d]];
+			s0 ^= gfmul[d0][gfcauchy[3][d]];
+		}
+
+		/* first disk with all coefficients at 1 */
+		d0 = v_8(v[0][i]);
+
+		p0 ^= d0;
+		q0 ^= d0;
+		r0 ^= d0;
+		s0 ^= d0;
+
+		v_8(p[i]) = p0;
+		v_8(q[i]) = q0;
+		v_8(r[i]) = r0;
+		v_8(s[i]) = s0;
+	}
+}
+
+/*
+ * PAR5 (penta parity with Cauchy matrix) 8bit C implementation
+ * 
+ * Note that instead of a generic multiplicationt table, likely resulting
+ * in multiple cache misses, a precomputed table could be used. 
+ * But this is only a kind of reference function, and we are not really
+ * interested in speed.
+ */
+void raid_par5_int8(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	uint8_t* r;
+	uint8_t* s;
+	uint8_t* t;
+	int d, l;
+	size_t i;
+
+	uint8_t d0, t0, s0, r0, q0, p0;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+	r = v[nd+2];
+	s = v[nd+3];
+	t = v[nd+4];
+
+	for(i=0;i<size;i+=1) {
+		p0 = q0 = r0 = s0 = t0 = 0;
+		for(d=l;d>0;--d) {
+			d0 = v_8(v[d][i]);
+
+			p0 ^= d0;
+			q0 ^= gfmul[d0][gfcauchy[1][d]];
+			r0 ^= gfmul[d0][gfcauchy[2][d]];
+			s0 ^= gfmul[d0][gfcauchy[3][d]];
+			t0 ^= gfmul[d0][gfcauchy[4][d]];
+		}
+
+		/* first disk with all coefficients at 1 */
+		d0 = v_8(v[0][i]);
+
+		p0 ^= d0;
+		q0 ^= d0;
+		r0 ^= d0;
+		s0 ^= d0;
+		t0 ^= d0;
+
+		v_8(p[i]) = p0;
+		v_8(q[i]) = q0;
+		v_8(r[i]) = r0;
+		v_8(s[i]) = s0;
+		v_8(t[i]) = t0;
+	}
+}
+
+/*
+ * PAR6 (hexa parity with Cauchy matrix) 8bit C implementation
+ * 
+ * Note that instead of a generic multiplicationt table, likely resulting
+ * in multiple cache misses, a precomputed table could be used. 
+ * But this is only a kind of reference function, and we are not really
+ * interested in speed.
+ */
+void raid_par6_int8(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	uint8_t* r;
+	uint8_t* s;
+	uint8_t* t;
+	uint8_t* u;
+	int d, l;
+	size_t i;
+
+	uint8_t d0, u0, t0, s0, r0, q0, p0;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+	r = v[nd+2];
+	s = v[nd+3];
+	t = v[nd+4];
+	u = v[nd+5];
+
+	for(i=0;i<size;i+=1) {
+		p0 = q0 = r0 = s0 = t0 = u0 = 0;
+		for(d=l;d>0;--d) {
+			d0 = v_8(v[d][i]);
+
+			p0 ^= d0;
+			q0 ^= gfmul[d0][gfcauchy[1][d]];
+			r0 ^= gfmul[d0][gfcauchy[2][d]];
+			s0 ^= gfmul[d0][gfcauchy[3][d]];
+			t0 ^= gfmul[d0][gfcauchy[4][d]];
+			u0 ^= gfmul[d0][gfcauchy[5][d]];
+		}
+
+		/* first disk with all coefficients at 1 */
+		d0 = v_8(v[0][i]);
+
+		p0 ^= d0;
+		q0 ^= d0;
+		r0 ^= d0;
+		s0 ^= d0;
+		t0 ^= d0;
+		u0 ^= d0;
+
+		v_8(p[i]) = p0;
+		v_8(q[i]) = q0;
+		v_8(r[i]) = r0;
+		v_8(s[i]) = s0;
+		v_8(t[i]) = t0;
+		v_8(u[i]) = u0;
+	}
+}
+
+/*
+ * Recover failure of one data block at index id[0] using parity at index
+ * ip[0] for any RAID level.
+ *
+ * Starting from the equation:
+ *
+ * Pd = A[ip[0],id[0]] * Dx
+ *
+ * and solving we get:
+ *
+ * Dx = A[ip[0],id[0]]^-1 * Pd
+ */
+void raid_rec1_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* pa;
+	const uint8_t* T;
+	uint8_t G;
+	uint8_t V;
+	size_t i;
+
+	(void)nr; /* unused, it's always 1 */
+
+	/* if it's RAID5 uses the dedicated and faster function */
+	if (ip[0] == 0) {
+		raid_rec1_par1(id, nd, size, vv);
+		return;
+	}
+
+	/* setup the coefficients matrix */
+	G = A(ip[0],id[0]);
+
+	/* invert it to solve the system of linear equations */
+	V = inv(G);
+
+	/* get multiplication tables */
+	T = table(V);
+
+	/* compute delta parity */
+	raid_delta_gen(1, id, ip, nd, size, vv);
+
+	p = v[nd+ip[0]];
+	pa = v[id[0]];
+
+	for(i=0;i<size;++i) {
+		/* delta */
+		uint8_t Pd = p[i] ^ pa[i];
+
+		/* reconstruct */
+		pa[i] = T[Pd];
+	}
+}
+
+/*
+ * Recover failure of two data blocks at indexes id[0],id[1] using parity at
+ * indexes ip[0],ip[1] for any RAID level.
+ *
+ * Starting from the equations:
+ *
+ * Pd = A[ip[0],id[0]] * Dx + A[ip[0],id[1]] * Dy
+ * Qd = A[ip[1],id[0]] * Dx + A[ip[1],id[1]] * Dy
+ *
+ * we solve inverting the coefficients matrix.
+ */
+void raid_rec2_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* pa;
+	uint8_t* q;
+	uint8_t* qa;
+	const int N = 2;
+	const uint8_t* T[N][N];
+	uint8_t G[N*N];
+	uint8_t V[N*N];
+	size_t i;
+	int j, k;
+
+	(void)nr; /* unused, it's always 2 */
+
+	/* if it's RAID6 uses the dedicated and faster function */
+	if (ip[0] == 0 && ip[1] == 1) {
+		raid_rec2_par2(id, ip, nd, size, vv);
+		return;
+	}
+
+	/* setup the coefficients matrix */
+	for(j=0;j<N;++j) {
+		for(k=0;k<N;++k) {
+			G[j*N+k] = A(ip[j],id[k]);
+		}
+	}
+
+	/* invert it to solve the system of linear equations */
+	raid_invert(G, V, N);
+
+	/* get multiplication tables */
+	for(j=0;j<N;++j) {
+		for(k=0;k<N;++k) {
+			T[j][k] = table( V[j*N+k] );
+		}
+	}
+
+	/* compute delta parity */
+	raid_delta_gen(2, id, ip, nd, size, vv);
+
+	p = v[nd+ip[0]];
+	q = v[nd+ip[1]];
+	pa = v[id[0]];
+	qa = v[id[1]];
+
+	for(i=0;i<size;++i) {
+		/* delta */
+		uint8_t Pd = p[i] ^ pa[i];
+		uint8_t Qd = q[i] ^ qa[i];
+
+		/* reconstruct */
+		pa[i] = T[0][0][Pd] ^ T[0][1][Qd];
+		qa[i] = T[1][0][Pd] ^ T[1][1][Qd];
+	}
+}
+
+/*
+ * Recover failure of N data blocks at indexes id[N] using parity at indexes
+ * ip[N] for any RAID level.
+ *
+ * Starting from the N equations, with 0<=i<N :
+ *
+ * PD[i] = sum(A[ip[i],id[j]] * D[i]) 0<=j<N
+ *
+ * we solve inverting the coefficients matrix.
+ *
+ * Note that referring at previous equations you have:
+ * PD[0] = Pd, PD[1] = Qd, PD[2] = Rd, ...
+ * D[0] = Dx, D[1] = Dy, D[2] = Dz, ...
+ */
+void raid_recX_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p[RAID_PARITY_MAX];
+	uint8_t* pa[RAID_PARITY_MAX];
+	const uint8_t* T[RAID_PARITY_MAX][RAID_PARITY_MAX];
+	uint8_t G[RAID_PARITY_MAX*RAID_PARITY_MAX];
+	uint8_t V[RAID_PARITY_MAX*RAID_PARITY_MAX];
+	size_t i;
+	int j, k;
+
+	/* setup the coefficients matrix */
+	for(j=0;j<nr;++j) {
+		for(k=0;k<nr;++k) {
+			G[j*nr+k] = A(ip[j],id[k]);
+		}
+	}
+
+	/* invert it to solve the system of linear equations */
+	raid_invert(G, V, nr);
+
+	/* get multiplication tables */
+	for(j=0;j<nr;++j) {
+		for(k=0;k<nr;++k) {
+			T[j][k] = table( V[j*nr+k] );
+		}
+	}
+
+	/* compute delta parity */
+	raid_delta_gen(nr, id, ip, nd, size, vv);
+
+	for(j=0;j<nr;++j) {
+		p[j] = v[nd+ip[j]];
+		pa[j] = v[id[j]];
+	}
+
+	for(i=0;i<size;++i) {
+		uint8_t PD[RAID_PARITY_MAX];
+
+		/* delta */
+		for(j=0;j<nr;++j)
+			PD[j] = p[j][i] ^ pa[j][i];
+
+		/* reconstruct */
+		for(j=0;j<nr;++j) {
+			uint8_t b = 0;
+			for(k=0;k<nr;++k)
+				b ^= T[j][k][PD[k]];
+			pa[j][i] = b;
+		}
+	}
+}
+
diff --git a/raid/internal.h b/raid/internal.h
new file mode 100644
index 00000000..e1ac2663
--- /dev/null
+++ b/raid/internal.h
@@ -0,0 +1,107 @@
+/*
+ * Copyright (C) 2013 Andrea Mazzoleni
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef __RAID_INTERNAL_H
+#define __RAID_INTERNAL_H
+
+#include "../portable.h"
+
+#include "raid.h"
+
+/*
+ * Internal functions.
+ *
+ * These are intented to provide access for testing.
+ */
+void raid_invert(uint8_t* M, uint8_t* V, int n);
+void raid_delta_gen(int nr, const int* id, const int* ip, int nd, size_t size, void** v);
+void raid_rec1_par1(const int* id, int nd, size_t size, void** v);
+void raid_rec2_par2(const int* id, const int* ip, int nd, size_t size, void** vv);
+void raid_par1_int32(int nd, size_t size, void** vv);
+void raid_par1_int64(int nd, size_t size, void** vv);
+void raid_par1_sse2(int nd, size_t size, void** vv);
+void raid_par2_int32(int nd, size_t size, void** vv);
+void raid_par2_int64(int nd, size_t size, void** vv);
+void raid_par2_sse2(int nd, size_t size, void** vv);
+void raid_par2_sse2ext(int nd, size_t size, void** vv);
+void raid_parz_int32(int nd, size_t size, void** vv);
+void raid_parz_int64(int nd, size_t size, void** vv);
+void raid_parz_sse2(int nd, size_t size, void** vv);
+void raid_parz_sse2ext(int nd, size_t size, void** vv);
+void raid_par3_int8(int nd, size_t size, void** vv);
+void raid_par3_ssse3(int nd, size_t size, void** vv);
+void raid_par3_ssse3ext(int nd, size_t size, void** vv);
+void raid_par4_int8(int nd, size_t size, void** vv);
+void raid_par4_ssse3(int nd, size_t size, void** vv);
+void raid_par4_ssse3ext(int nd, size_t size, void** vv);
+void raid_par5_int8(int nd, size_t size, void** vv);
+void raid_par5_ssse3(int nd, size_t size, void** vv);
+void raid_par5_ssse3ext(int nd, size_t size, void** vv);
+void raid_par6_int8(int nd, size_t size, void** vv);
+void raid_par6_ssse3(int nd, size_t size, void** vv);
+void raid_par6_ssse3ext(int nd, size_t size, void** vv);
+void raid_rec1_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
+void raid_rec2_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
+void raid_recX_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
+void raid_rec1_ssse3(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
+void raid_rec2_ssse3(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
+void raid_recX_ssse3(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
+
+/*
+ * Internal naming.
+ *
+ * These are intented to provide access for testing.
+ */
+const char* raid_par1_tag(void);
+const char* raid_par2_tag(void);
+const char* raid_parz_tag(void);
+const char* raid_par3_tag(void);
+const char* raid_par4_tag(void);
+const char* raid_par5_tag(void);
+const char* raid_par6_tag(void);
+const char* raid_rec1_tag(void);
+const char* raid_rec2_tag(void);
+const char* raid_recX_tag(void);
+
+/*
+ * Internal forwarders.
+ */
+extern void (*raid_par_ptr[RAID_PARITY_MAX])(int nd, size_t size, void** vv);
+extern void (*raid_par3_ptr)(int nd, size_t size, void** vv);
+extern void (*raid_parz_ptr)(int nd, size_t size, void** vv); 
+extern void (*raid_rec_ptr[RAID_PARITY_MAX])(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
+
+/*
+ * Tables.
+ */
+extern const uint8_t raid_gfmul[256][256] __attribute__((aligned(256)));
+extern const uint8_t raid_gfexp[256] __attribute__((aligned(256)));
+extern const uint8_t raid_gfinv[256] __attribute__((aligned(256)));
+extern const uint8_t raid_gfvandermonde[3][256] __attribute__((aligned(256)));
+extern const uint8_t raid_gfcauchy[6][256] __attribute__((aligned(256)));
+extern const uint8_t raid_gfcauchypshufb[251][4][2][16] __attribute__((aligned(256)));
+extern const uint8_t raid_gfmulpshufb[256][2][16] __attribute__((aligned(256)));
+#define gfmul raid_gfmul
+#define gfexp raid_gfexp
+#define gfinv raid_gfinv
+#define gfvandermonde raid_gfvandermonde
+#define gfcauchy raid_gfcauchy
+#define gfcauchypshufb raid_gfcauchypshufb
+#define gfmulpshufb raid_gfmulpshufb
+
+#endif
+
diff --git a/raid/intz.c b/raid/intz.c
new file mode 100644
index 00000000..f412f001
--- /dev/null
+++ b/raid/intz.c
@@ -0,0 +1,122 @@
+/*
+ * Copyright (C) 2013 Andrea Mazzoleni
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "internal.h"
+#include "gf.h"
+
+/*
+ * PARz (triple parity with powers of 2^-1) 32bit C implementation
+ */
+void raid_parz_int32(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	uint8_t* r;
+	int d, l;
+	size_t i;
+
+	uint32_t d0, r0, q0, p0;
+	uint32_t d1, r1, q1, p1;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+	r = v[nd+2];
+
+	for(i=0;i<size;i+=8) {
+		r0 = q0 = p0 = v_32(v[l][i]);
+		r1 = q1 = p1 = v_32(v[l][i+4]);
+		for(d=l-1;d>=0;--d) {
+			d0 = v_32(v[d][i]);
+			d1 = v_32(v[d][i+4]);
+
+			p0 ^= d0;
+			p1 ^= d1;
+
+			q0 = x2_32(q0);
+			q1 = x2_32(q1);
+
+			q0 ^= d0;
+			q1 ^= d1;
+
+			r0 = d2_32(r0);
+			r1 = d2_32(r1);
+
+			r0 ^= d0;
+			r1 ^= d1;
+		}
+		v_32(p[i]) = p0;
+		v_32(p[i+4]) = p1;
+		v_32(q[i]) = q0;
+		v_32(q[i+4]) = q1;
+		v_32(r[i]) = r0;
+		v_32(r[i+4]) = r1;
+	}
+}
+
+/*
+ * PARz (triple parity with powers of 2^-1) 64bit C implementation
+ */
+void raid_parz_int64(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	uint8_t* r;
+	int d, l;
+	size_t i;
+
+	uint64_t d0, r0, q0, p0;
+	uint64_t d1, r1, q1, p1;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+	r = v[nd+2];
+
+	for(i=0;i<size;i+=16) {
+		r0 = q0 = p0 = v_64(v[l][i]);
+		r1 = q1 = p1 = v_64(v[l][i+8]);
+		for(d=l-1;d>=0;--d) {
+			d0 = v_64(v[d][i]);
+			d1 = v_64(v[d][i+8]);
+
+			p0 ^= d0;
+			p1 ^= d1;
+
+			q0 = x2_64(q0);
+			q1 = x2_64(q1);
+
+			q0 ^= d0;
+			q1 ^= d1;
+
+			r0 = d2_64(r0);
+			r1 = d2_64(r1);
+
+			r0 ^= d0;
+			r1 ^= d1;
+		}
+		v_64(p[i]) = p0;
+		v_64(p[i+8]) = p1;
+		v_64(q[i]) = q0;
+		v_64(q[i+8]) = q1;
+		v_64(r[i]) = r0;
+		v_64(r[i+8]) = r1;
+	}
+}
+
diff --git a/mktables.c b/raid/mktables.c
similarity index 92%
rename from mktables.c
rename to raid/mktables.c
index 2f726d5d..932ba1fb 100644
--- a/mktables.c
+++ b/raid/mktables.c
@@ -201,8 +201,11 @@ int main(void)
 	printf(" */\n");
 	printf("\n");
 
+	printf("#include \"internal.h\"\n");
+	printf("\n");
+
 	/* a*b */
-	printf("const uint8_t  __attribute__((aligned(256))) gfmul[256][256] =\n");
+	printf("const uint8_t __attribute__((aligned(256))) raid_gfmul[256][256] =\n");
 	printf("{\n");
 	for(i=0;i<256;++i) {
 		printf("\t{\n");
@@ -223,7 +226,7 @@ int main(void)
 	printf("};\n\n");
 
 	/* 2^a */
-	printf("const uint8_t __attribute__((aligned(256))) gfexp[256] =\n");
+	printf("const uint8_t __attribute__((aligned(256))) raid_gfexp[256] =\n");
 	printf("{\n");
 	v = 1;
 	for(i=0;i<256;++i) {
@@ -239,7 +242,7 @@ int main(void)
 	printf("};\n\n");
 
 	/* 1/a */
-	printf("const uint8_t __attribute__((aligned(256))) gfinv[256] =\n");
+	printf("const uint8_t __attribute__((aligned(256))) raid_gfinv[256] =\n");
 	printf("{\n");
 	printf("\t/* note that the first element is not significative */\n");
 	for(i=0;i<256;++i) {
@@ -272,7 +275,7 @@ int main(void)
 		printf("\n");
 	}
 	printf(" */\n");
-	printf("const uint8_t  __attribute__((aligned(256))) gfvandermonde[%u][256] =\n", 3);
+	printf("const uint8_t __attribute__((aligned(256))) raid_gfvandermonde[%u][256] =\n", 3);
 	printf("{\n");
 	for(p=0;p<3;++p) {
 		printf("\t{\n");
@@ -304,7 +307,7 @@ int main(void)
 		printf("\n");
 	}
 	printf(" */\n");
-	printf("const uint8_t  __attribute__((aligned(256))) gfcauchy[%u][256] =\n", PARITY);
+	printf("const uint8_t __attribute__((aligned(256))) raid_gfcauchy[%u][256] =\n", PARITY);
 	printf("{\n");
 	for(p=0;p<PARITY;++p) {
 		printf("\t{\n");
@@ -328,7 +331,7 @@ int main(void)
 	printf(" * Indexes are [DISK][PARITY - 2][LH].\n");
 	printf(" * Where DISK is from 0 to %u, PARITY from 2 to %u, LH from 0 to 1.\n", DISK - 1, PARITY - 1);
 	printf(" */\n");
-	printf("const uint8_t  __attribute__((aligned(256))) gfcauchypshufb[%u][%u][2][16] =\n", DISK, np(PARITY - 2));
+	printf("const uint8_t __attribute__((aligned(256))) raid_gfcauchypshufb[%u][%u][2][16] =\n", DISK, np(PARITY - 2));
 	printf("{\n");
 	for(i=0;i<DISK;++i) {
 		printf("\t{\n");
@@ -360,7 +363,7 @@ int main(void)
 	printf(" * Indexes are [MULTIPLER][LH].\n");
 	printf(" * Where MULTIPLER is from 0 to 255, LH from 0 to 1.\n");
 	printf(" */\n");
-	printf("const uint8_t  __attribute__((aligned(256))) gfmulpshufb[256][2][16] =\n");
+	printf("const uint8_t __attribute__((aligned(256))) raid_gfmulpshufb[256][2][16] =\n");
 	printf("{\n");
 	for(i=0;i<256;++i) {
 		printf("\t{\n");
diff --git a/raid/module.c b/raid/module.c
new file mode 100644
index 00000000..db83bd16
--- /dev/null
+++ b/raid/module.c
@@ -0,0 +1,189 @@
+/*
+ * Copyright (C) 2013 Andrea Mazzoleni
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "internal.h"
+#include "cpu.h"
+
+void raid_init(void)
+{
+	raid_par3_ptr = raid_par3_int8;
+	raid_par_ptr[3] = raid_par4_int8;
+	raid_par_ptr[4] = raid_par5_int8;
+	raid_par_ptr[5] = raid_par6_int8;
+
+	if (sizeof(void*) == 4) {
+		raid_par_ptr[0] = raid_par1_int32;
+		raid_par_ptr[1] = raid_par2_int32;
+		raid_parz_ptr = raid_parz_int32;
+	} else {
+		raid_par_ptr[0] = raid_par1_int64;
+		raid_par_ptr[1] = raid_par2_int64;
+		raid_parz_ptr = raid_parz_int64;
+	}
+
+	raid_rec_ptr[0] = raid_rec1_int8;
+	raid_rec_ptr[1] = raid_rec2_int8;
+	raid_rec_ptr[2] = raid_recX_int8;
+	raid_rec_ptr[3] = raid_recX_int8;
+	raid_rec_ptr[4] = raid_recX_int8;
+	raid_rec_ptr[5] = raid_recX_int8;
+
+#if defined(__i386__) || defined(__x86_64__)
+	if (raid_cpu_has_sse2()) {
+		raid_par_ptr[0] = raid_par1_sse2;
+#if defined(__x86_64__)
+		if (raid_cpu_has_slowextendedreg()) {
+			raid_par_ptr[1] = raid_par2_sse2;
+			raid_parz_ptr = raid_parz_sse2;
+		} else {
+			raid_par_ptr[1] = raid_par2_sse2ext;
+			raid_parz_ptr = raid_parz_sse2ext;
+		}
+#else
+		raid_par_ptr[1] = raid_par2_sse2;
+		raid_parz_ptr = raid_parz_sse2;
+#endif
+	}
+
+	if (raid_cpu_has_ssse3()) {
+#if defined(__x86_64__)
+		raid_par3_ptr = raid_par3_ssse3ext;
+		raid_par_ptr[3] = raid_par4_ssse3ext;
+		raid_par_ptr[4] = raid_par5_ssse3ext;
+		raid_par_ptr[5] = raid_par6_ssse3ext;
+#else
+		raid_par3_ptr = raid_par3_ssse3;
+		raid_par_ptr[3] = raid_par4_ssse3;
+		raid_par_ptr[4] = raid_par5_ssse3;
+		raid_par_ptr[5] = raid_par6_ssse3;
+#endif
+		raid_rec_ptr[0] = raid_rec1_ssse3;
+		raid_rec_ptr[1] = raid_rec2_ssse3;
+		raid_rec_ptr[2] = raid_recX_ssse3;
+		raid_rec_ptr[3] = raid_recX_ssse3;
+		raid_rec_ptr[4] = raid_recX_ssse3;
+		raid_rec_ptr[5] = raid_recX_ssse3;
+	}
+#endif
+
+	/* set the default mode */
+	raid_mode(RAID_MODE_CAUCHY);
+}
+
+static struct raid_func {
+	const char* name;
+	void* p;
+} RAID_FUNC[] = {
+	{ "int8", raid_par3_int8 },
+	{ "int8", raid_par4_int8 },
+	{ "int8", raid_par5_int8 },
+	{ "int8", raid_par6_int8 },
+	{ "int32", raid_par1_int32 },
+	{ "int64", raid_par1_int64 },
+	{ "int32", raid_par2_int32 },
+	{ "int64", raid_par2_int64 },
+	{ "int32", raid_parz_int32 },
+	{ "int64", raid_parz_int64 },
+	{ "int8", raid_rec1_int8 },
+	{ "int8", raid_rec2_int8 },
+	{ "int8", raid_recX_int8 },
+
+#if defined(__i386__) || defined(__x86_64__)
+	{ "sse2", raid_par1_sse2 },
+	{ "sse2", raid_par2_sse2 },
+	{ "sse2", raid_parz_sse2 },
+	{ "ssse3", raid_par3_ssse3 },
+	{ "ssse3", raid_par4_ssse3 },
+	{ "ssse3", raid_par5_ssse3 },
+	{ "ssse3", raid_par6_ssse3 },
+	{ "ssse3", raid_rec1_ssse3 },
+	{ "ssse3", raid_rec2_ssse3 },
+	{ "ssse3", raid_recX_ssse3 },
+#endif
+
+#if defined(__x86_64__)
+	{ "sse2e", raid_par2_sse2ext },
+	{ "sse2e", raid_parz_sse2ext },
+	{ "ssse3e", raid_par3_ssse3ext },
+	{ "ssse3e", raid_par4_ssse3ext },
+	{ "ssse3e", raid_par5_ssse3ext },
+	{ "ssse3e", raid_par6_ssse3ext },
+#endif
+	{ 0, 0 }
+};
+
+static const char* raid_tag(void* func)
+{
+	struct raid_func* i = RAID_FUNC;
+	while (i->name != 0) {
+		if (i->p == func)
+			return i->name;
+		++i;
+	}
+	return "unknown";
+}
+
+const char* raid_par1_tag(void)
+{
+	return raid_tag(raid_par_ptr[0]);
+}
+
+const char* raid_par2_tag(void)
+{
+	return raid_tag(raid_par_ptr[1]);
+}
+
+const char* raid_parz_tag(void)
+{
+	return raid_tag(raid_parz_ptr);
+}
+
+const char* raid_par3_tag(void)
+{
+	return raid_tag(raid_par_ptr[2]);
+}
+
+const char* raid_par4_tag(void)
+{
+	return raid_tag(raid_par_ptr[3]);
+}
+
+const char* raid_par5_tag(void)
+{
+	return raid_tag(raid_par_ptr[4]);
+}
+
+const char* raid_par6_tag(void)
+{
+	return raid_tag(raid_par_ptr[5]);
+}
+
+const char* raid_rec1_tag(void)
+{
+	return raid_tag(raid_rec_ptr[0]);
+}
+
+const char* raid_rec2_tag(void)
+{
+	return raid_tag(raid_rec_ptr[1]);
+}
+
+const char* raid_recX_tag(void)
+{
+	return raid_tag(raid_rec_ptr[2]);
+}
+
diff --git a/raid/raid.c b/raid/raid.c
new file mode 100644
index 00000000..bdad2a26
--- /dev/null
+++ b/raid/raid.c
@@ -0,0 +1,458 @@
+/*
+ * Copyright (C) 2013 Andrea Mazzoleni
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/*
+ * Copyright notes:
+ *
+ * The RAID5 and RAID6 support was originally derived from the
+ * H. Peter Anvin paper "The mathematics of RAID-6" [1] and the
+ * libraid6 library, also by H. Peter Anvin, released with license
+ * "GPLv2 or any later version" inside the Linux Kernel 2.6.38.
+ *
+ * This support was later completely rewritten (many times), but
+ * the H. Peter Anvin's Copyright may still apply.
+ *
+ * The others RAID levels and the recovering based on matrix
+ * inversion is original work implemented from scratch.
+ */
+
+/*
+ * The RAID5 and RAID6 support is implemented using the Galois Field
+ * GF(2^8) with the primitive polynomial x^8 + x^4 + x^3 + x^2 + 1 
+ * (285 decimal).
+ *
+ * The parity P and Q of a set of N disk Di with 0<=i<N, is computed
+ * using the equations:
+ *
+ * P = sum(Di)
+ * Q = sum(2^i * Di) with 0<=i<N
+ *
+ * This approach is the same used by the Linux Kernel RAID, and by ZFS
+ * RAIDZ2, better described in the H. Peter Anvin paper "The mathematics 
+ * of RAID-6" [1]. 
+ *
+ * To support triple parity, it was first evaluated and then dropped, an
+ * extension of the same approach, with additional parity coefficients set
+ * as powers of 2^-1, with equations:
+ *
+ * P = sum(Di)
+ * Q = sum(2^i * Di)
+ * R = sum(2^-i * Di) with 0<=i<N
+ *
+ * This approach works well for triple parity and it's very efficient,
+ * because we can implement very fast parallel multiplications and
+ * divisions by 2 in GF(2^8).
+ *
+ * It's also similar at the approach used by ZFS RAIDZ3, with the 
+ * difference that ZFS uses powers of 4 instead of 2^-1.
+ *
+ * Unfortunately it doesn't work beyond triple parity, because whatever
+ * value we choose to generate the power coefficients to compute other
+ * parities, the resulting equations are not solvable for some
+ * combinations of missing disks.
+ *
+ * This is expected, because the Vandermonde matrix used to compute the
+ * parity has no guarantee to have all submatrices not singular
+ * [2, Chap 11, Problem 7] and this is a requirement to have
+ * a MDS (Maximum Distance Separable) code [2, Chap 11, Theorem 8].
+ *
+ * To overcome this limitation, we use a Cauchy matrix [3][4] to compute 
+ * the parity. A Cauchy matrix has the property to have all the square
+ * submatrices not singular, resulting in always solvable equations,
+ * for any combination of missing disks.
+ *
+ * The problem of this approach is that it requires the use of
+ * generic multiplications, and not only by 2 or 2^-1, potentially 
+ * affecting badly the performance.
+ *
+ * Hopefully there is a method to implement parallel multiplications
+ * using SSSE3 instructions [1][5]. Method competitive with the
+ * computation of triple parity using power coefficients.
+ *
+ * Another important property of the Cauchy matrix is that we can setup
+ * the first two rows with coeffients equal at the RAID5 and RAID6 approach
+ * decribed, resulting in a compatible extension, and requiring SSSE3
+ * instructions only if triple parity or beyond is used.
+ * 
+ * The matrix is also adjusted, multipling each row by a constant factor
+ * to make the first column of all 1, to optimize the computation for
+ * the first disk.
+ *
+ * This results in the matrix A[row,col] defined as:
+ * 
+ * 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01...
+ * 01 02 04 08 10 20 40 80 1d 3a 74 e8 cd 87 13 26 4c 98 2d 5a b4 75...
+ * 01 f5 d2 c4 9a 71 f1 7f fc 87 c1 c6 19 2f 40 55 3d ba 53 04 9c 61...
+ * 01 bb a6 d7 c7 07 ce 82 4a 2f a5 9b b6 60 f1 ad e7 f4 06 d2 df 2e...
+ * 01 97 7f 9c 7c 18 bd a2 58 1a da 74 70 a3 e5 47 29 07 f5 80 23 e9...
+ * 01 2b 3f cf 73 2c d6 ed cb 74 15 78 8a c1 17 c9 89 68 21 ab 76 3b...
+ * (see tables.h for the matrix with all the columns)
+ *
+ * This matrix supports 6 level of parity, one for each row, for up to 251
+ * data disks, one for each column, with all the 377,342,351,231 square
+ * submatrices not singular, verified also with brute-force.
+ * 
+ * This matrix can be extended to support any number of parities, just 
+ * adding additional rows, and removing one column for each new row.
+ * (see mktables.c for more details in how the matrix is generated)
+ * 
+ * In details, parity is computed as:
+ *
+ * P = sum(Di)
+ * Q = sum(2^i *  Di)
+ * R = sum(A[2,i] * Di)
+ * S = sum(A[3,i] * Di)
+ * T = sum(A[4,i] * Di)
+ * U = sum(A[5,i] * Di) with 0<=i<N
+ *
+ * To recover from a failure of six disks at indexes x,y,z,h,v,w,
+ * with 0<=x<y<z<h<v<w<N, we compute the parity of the available N-6 
+ * disks as:
+ *
+ * Pa = sum(Di)
+ * Qa = sum(2^i * Di)
+ * Ra = sum(A[2,i] * Di)
+ * Sa = sum(A[3,i] * Di)
+ * Ta = sum(A[4,i] * Di)
+ * Ua = sum(A[5,i] * Di) with 0<=i<N,i!=x,i!=y,i!=z,i!=h,i!=v,i!=w.
+ *
+ * And if we define:
+ *
+ * Pd = Pa + P
+ * Qd = Qa + Q
+ * Rd = Ra + R
+ * Sd = Sa + S
+ * Td = Ta + T
+ * Ud = Ua + U
+ *
+ * we can sum these two sets of equations, obtaining:
+ *
+ * Pd =          Dx +          Dy +          Dz +          Dh +          Dv +          Dw
+ * Qd =    2^x * Dx +    2^y * Dy +    2^z * Dz +    2^h * Dh +    2^v * Dv +    2^w * Dw
+ * Rd = A[2,x] * Dx + A[2,y] * Dy + A[2,z] * Dz + A[2,h] * Dh + A[2,v] * Dv + A[2,w] * Dw
+ * Sd = A[3,x] * Dx + A[3,y] * Dy + A[3,z] * Dz + A[3,h] * Dh + A[3,v] * Dv + A[3,w] * Dw
+ * Td = A[4,x] * Dx + A[4,y] * Dy + A[4,z] * Dz + A[4,h] * Dh + A[4,v] * Dv + A[4,w] * Dw
+ * Ud = A[5,x] * Dx + A[5,y] * Dy + A[5,z] * Dz + A[5,h] * Dh + A[5,v] * Dv + A[5,w] * Dw
+ *
+ * A linear system always solvable because the coefficients matrix is
+ * always not singular due the properties of the matrix A[].
+ *
+ * Resulting speed in x64, with 8 data disks, using a stripe of 256 KiB,
+ * for a Core i7-3740QM CPU @ 2.7GHz is:
+ * 
+ *             int8   int32   int64    sse2   sse2e   ssse3  ssse3e
+ *   par1             11927   22075   36004
+ *   par2              3378    5874   18235   19164
+ *   par3       844                                    8814    9419
+ *   par4       665                                    6836    7415
+ *   par5       537                                    5388    5686
+ *   par6       449                                    4307    4789
+ *
+ * Values are in MiB/s of data processed, not counting generated parity.
+ *
+ * You can replicate these results in your machine using the "snapraid -T"
+ * command.
+ *
+ * For comparison, the triple parity computation using the power
+ * coeffients "1,2,2^-1" is only a little faster than the one based on
+ * the Cauchy matrix if SSSE3 is present.
+ *
+ *             int8   int32   int64    sse2   sse2e   ssse3  ssse3e
+ *   parz              2112    3118    9589   10304
+ *
+ *
+ * In conclusion, the use of power coefficients, and specifically powers
+ * of 1,2,2^-1, is the best option to implement triple parity in CPUs 
+ * without SSSE3.
+ * But if a modern CPU with SSSE3 (or similar) is available, the Cauchy
+ * matrix is the best option because it provides a fast and general 
+ * approach working for any number of parities.
+ *
+ * References:
+ * [1] Anvin, "The mathematics of RAID-6", 2004
+ * [2] MacWilliams, Sloane, "The Theory of Error-Correcting Codes", 1977
+ * [3] Blomer, "An XOR-Based Erasure-Resilient Coding Scheme", 1995
+ * [4] Roth, "Introduction to Coding Theory", 2006
+ * [5] Plank, "Screaming Fast Galois Field Arithmetic Using Intel SIMD Instructions", 2013
+ */
+
+#include "internal.h"
+#include "gf.h"
+
+/**
+ * Generator matrix currently used.
+ */
+const uint8_t (*raid_gen)[256];
+
+void raid_mode(int mode)
+{
+	if (mode == RAID_MODE_VANDERMONDE) {
+		raid_par_ptr[2] = raid_parz_ptr;
+		raid_gen = gfvandermonde;
+	} else {
+		raid_par_ptr[2] = raid_par3_ptr;
+		raid_gen = gfcauchy;
+	}
+}
+
+/**
+ * Buffer filled with 0 used in recovering.
+ */
+static void* raid_empty_zero_block;
+
+void raid_zero(void* zero)
+{
+	raid_empty_zero_block = zero;
+}
+
+/* internal forwarder */
+void (*raid_par_ptr[RAID_PARITY_MAX])(int nd, size_t size, void** vv);
+void (*raid_par3_ptr)(int nd, size_t size, void** vv);
+void (*raid_parz_ptr)(int nd, size_t size, void** vv);
+
+void raid_par(int np, int nd, size_t size, void** v)
+{
+	BUG_ON(np < 1 || np > RAID_PARITY_MAX);
+
+	raid_par_ptr[np - 1](nd, size, v);
+}
+
+/**
+ * Inverts the square matrix M of size nxn into V.
+ * We use Gauss elimination to invert.
+ */
+void raid_invert(uint8_t* M, uint8_t* V, int n)
+{
+	int i,j,k;
+
+	/* set the identity matrix in V */
+	for(i=0;i<n;++i) {
+		for(j=0;j<n;++j) {
+			V[i*n+j] = i == j;
+		}
+	}
+
+	/* for each element in the diagonal */
+	for(k=0;k<n;++k) {
+		uint8_t f;
+
+		/* the diagonal element cannot be 0 because */
+		/* we are inverting matrices with all the square submatrices */
+		/* not singular */
+		BUG_ON(M[k*n+k] == 0);
+
+		/* make the diagonal element to be 1 */
+		f = inv(M[k*n+k]);
+		for(j=0;j<n;++j) {
+			M[k*n+j] = mul(f, M[k*n+j]);
+			V[k*n+j] = mul(f, V[k*n+j]);
+		}
+
+		/* make all the elements over and under the diagonal to be 0 */
+		for(i=0;i<n;++i) {
+			if (i == k)
+				continue;
+			f = M[i*n+k];
+			for(j=0;j<n;++j) {
+				M[i*n+j] ^= mul(f, M[k*n+j]);
+				V[i*n+j] ^= mul(f, V[k*n+j]);
+			}
+		}
+	}
+}
+
+/**
+ * Computes the parity without the missing data blocks
+ * and store it in the buffers of such data blocks.
+ *
+ * This is the parity expressed as Pa,Qa,Ra,Sa,Ta,Ua
+ * in the equations.
+ *
+ * Note that all the other parities not in the ip[] vector
+ * are destroyed.
+ */
+void raid_delta_gen(int nr, const int* id, const int* ip, int nd, size_t size, void** v)
+{
+	void* p[RAID_PARITY_MAX];
+	void* pa[RAID_PARITY_MAX];
+	int i;
+
+	for(i=0;i<nr;++i) {
+		/* keep a copy of the parity buffer */
+		p[i] = v[nd+ip[i]];
+
+		/* buffer for missing data blocks */
+		pa[i] = v[id[i]];
+
+		/* set at zero the missing data blocks */
+		v[id[i]] = raid_empty_zero_block;
+
+		/* compute the parity over the missing data blocks */
+		v[nd+ip[i]] = pa[i];
+	}
+
+	/* recompute the minimal parity required */
+	raid_par(ip[nr - 1] + 1, nd, size, v);
+
+	for(i=0;i<nr;++i) {
+		/* restore disk buffers as before */
+		v[id[i]] = pa[i];
+
+		/* restore parity buffers as before */
+		v[nd+ip[i]] = p[i];
+	}
+}
+
+/**
+ * Recover failure of one data block for PAR1.
+ *
+ * Starting from the equation:
+ *
+ * Pd = Dx
+ *
+ * and solving we get:
+ *
+ * Dx = Pd
+ */
+void raid_rec1_par1(const int* id, int nd, size_t size, void** v)
+{
+	void* p;
+	void* pa;
+
+	/* for PAR1 we can directly compute the missing block */
+	/* and we don't need to use the zero buffer */
+	p = v[nd];
+	pa = v[id[0]];
+
+	/* use the parity as missing data block */
+	v[id[0]] = p;
+
+	/* compute the parity over the missing data block */
+	v[nd] = pa;
+
+	/* compute */
+	raid_par(1, nd, size, v);
+
+	/* restore as before */
+	v[id[0]] = pa;
+	v[nd] = p;
+}
+
+/**
+ * Recover failure of two data blocks for PAR2.
+ *
+ * Starting from the equations:
+ *
+ * Pd = Dx + Dy
+ * Qd = 2^id[0] * Dx + 2^id[1] * Dy
+ *
+ * and solving we get:
+ *
+ *               1                     2^(-id[0])
+ * Dy = ------------------- * Pd + ------------------- * Qd
+ *      2^(id[1]-id[0]) + 1        2^(id[1]-id[0]) + 1
+ *
+ * Dx = Dy + Pd
+ *
+ * with conditions:
+ *
+ * 2^id[0] != 0
+ * 2^(id[1]-id[0]) + 1 != 0
+ *
+ * That are always satisfied for any 0<=id[0]<id[1]<255.
+ */
+void raid_rec2_par2(const int* id, const int* ip, int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	size_t i;
+	uint8_t* p;
+	uint8_t* pa;
+	uint8_t* q;
+	uint8_t* qa;
+	const uint8_t* T[2];
+
+	/* get multiplication tables */
+	T[0] = table( inv(pow2(id[1]-id[0]) ^ 1) );
+	T[1] = table( inv(pow2(id[0]) ^ pow2(id[1])) );
+
+	/* compute delta parity */
+	raid_delta_gen(2, id, ip, nd, size, vv);
+
+	p = v[nd];
+	q = v[nd+1];
+	pa = v[id[0]];
+	qa = v[id[1]];
+
+	for(i=0;i<size;++i) {
+		/* delta */
+		uint8_t Pd = p[i] ^ pa[i];
+		uint8_t Qd = q[i] ^ qa[i];
+
+		/* reconstruct */
+		uint8_t Dy = T[0][Pd] ^ T[1][Qd];
+		uint8_t Dx = Pd ^ Dy;
+
+		/* set */
+		pa[i] = Dx;
+		qa[i] = Dy;
+	}
+}
+
+/* internal forwarder */
+void (*raid_rec_ptr[RAID_PARITY_MAX])(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
+
+void raid_rec(int nr, const int* id, const int* ip, int nd, size_t size, void** v)
+{
+	BUG_ON(nr > nd);
+
+	/* if failed data is present */
+	if (nr != 0) {
+		raid_rec_ptr[nr - 1](nr, id, ip, nd, size, v);
+	}
+}
+
+void raid_recpar(int nrd, const int* id, int nrp, int* ip, int np, int nd, size_t size, void** v)
+{
+	BUG_ON(nrd > nd);
+	BUG_ON(nrd + nrp > np);
+
+	/* if failed data is present */
+	if (nrd != 0) {
+		int iu[RAID_PARITY_MAX];
+		int i, j, k;
+
+		/* setup the vector of parities to use */
+		for(i=0,j=0,k=0;i<np;++i) {
+			if (ip[j] == i) {
+				++j;
+			} else {
+				iu[k] = i;
+				++k;
+			}
+		}
+
+		/* recover the data, and limit the parity use to needed ones */
+		raid_rec_ptr[nrd - 1](nrd, id, iu, nd, size, v);
+	}
+
+	/* recompute all the parities up to the last bad one */
+	if (nrp != 0) {
+		raid_par(ip[nrp - 1] + 1, nd, size, v);
+	}
+}
+
+
diff --git a/raid.h b/raid/raid.h
similarity index 55%
rename from raid.h
rename to raid/raid.h
index 78b7bee3..bbab02b4 100644
--- a/raid.h
+++ b/raid/raid.h
@@ -18,9 +18,6 @@
 #ifndef __RAID_H
 #define __RAID_H
 
-/****************************************************************************/
-/* raid */
-
 /**
  * RAID mode supporting up to 6 parities.
  *
@@ -62,8 +59,8 @@ void raid_init(void);
  * You can change mode at any time, and it will affect next calls to raid_par(),
  * raid_rec() and raid_recpar().
  *
- * The two modes are compatible for the first two levels of parity. The third one
- * is different.
+ * The two modes are compatible for the first two levels of parity.
+ * The third one is different.
  */
 void raid_mode(int mode);
 
@@ -80,9 +77,10 @@ void raid_zero(void* zero);
  *
  * \param np Number of parities to compute.
  * \param nd Number of data disks.
- * \param size Size of the blocks pointed by vv. It must be a power of 2, and at least 64.
+ * \param size Size of the blocks pointed by vv. It must be a multipler of 64.
  * \param vv Vector of pointers to the blocks for disks and parities.
- *   It has (::nd + ::np) elements. Each element points to a buffer of ::size bytes.
+ *   It has (::nd + ::np) elements.
+ *   Each element points to a buffer of ::size bytes.
  */
 void raid_par(int np, int nd, size_t size, void** vv);
 
@@ -93,12 +91,15 @@ void raid_par(int np, int nd, size_t size, void** vv);
  * If you want also to recompute parities, you can use raid_recpar().
  *
  * \param nr Number of failed data disks to recover.
- * \param id[] Vector of ::nr indexes of the data disks to recover. Starting from 0.
- * \param ip[] Vector of ::nr indexes of the parity disks to use for recovering. Starting from 0.
+ * \param id[] Vector of ::nr indexes of the data disks to recover.
+ *   The indexes start from 0. They must be in order.
+ * \param ip[] Vector of ::nr indexes of the parity disks to use for recovering.
+ *   The indexes start from 0. They must be in order.
  * \param nd Number of data disks.
- * \param size Size of the blocks pointed by vv. It must be a power of 2, and at least 64.
+ * \param size Size of the blocks pointed by vv. It must be a multipler of 64.
  * \param vv Vector of pointers to the blocks for disks and parities.
- *   It has (::nd + ::ip[::nr - 1] + 1) elements. Each element points to a buffer of ::size bytes.
+ *   It has (::nd + ::ip[::nr - 1] + 1) elements. The first elements are the blocks
+ *   for data, following with the parity blocks. Each blocks has ::size bytes.
  */
 void raid_rec(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
 
@@ -113,66 +114,21 @@ void raid_rec(int nr, const int* id, const int* ip, int nd, size_t size, void**
  * Ensure to have ::nrd + ::nrp <= ::np, otherwise recovering is not possible.
  *
  * \param nrd Number of failed data disks to recover.
- * \param id[] Vector of ::nrd indexes of the data disks to recover. Starting from 0.
+ * \param id[] Vector of ::nrd indexes of the data disks to recover.
+ *   The indexes start from 0. They must be in order.
  * \param nrp Number of failed parity disks to recover.
- * \param ip[] Vector of ::nrp indexes of the parity disks to recover. Starting from 0.
+ * \param ip[] Vector of ::nrp indexes of the parity disks to recover.
+ *   The indexes start from 0. They must be in order.
  *   Note that as difference than raid_rec() these are the failed parities.
- *   All the parities not specified here are assumed correct, and they are not recomputed.
+ *   All the parities not specified here are assumed correct, and they are
+ *   not recomputed.
  * \param np Number of parity disks.
  * \param nd Number of data disks.
- * \param size Size of the blocks pointed by vv. It must be a power of 2, and at least 64.
+ * \param size Size of the blocks pointed by vv. It must be a multipler of 64.
  * \param vv Vector of pointers to the blocks for disks and parities.
- *   It has (::nd + ::np) elements. Each element points to a buffer of ::size bytes.
+ *   It has (::nd + ::np) elements. The first elements are the blocks
+ *   for data, following with the parity blocks. Each blocks has ::size bytes. 
  */
 void raid_recpar(int nrd, const int* id, int nrp, int* ip, int np, int nd, size_t size, void** vv);
 
-/**
- * Gets the name of the selected function to compute the parity.
- * Only for testing.
- */
-const char* raid_par1_tag(void);
-const char* raid_par2_tag(void);
-const char* raid_parz_tag(void);
-const char* raid_par3_tag(void);
-const char* raid_par4_tag(void);
-const char* raid_par5_tag(void);
-const char* raid_par6_tag(void);
-const char* raid_rec1_tag(void);
-const char* raid_rec2_tag(void);
-const char* raid_recX_tag(void);
-
-/**
- * Specialized parity computation.
- * Only for testing.
- */
-void raid_par1_int32(int nd, size_t size, void** vv);
-void raid_par1_int64(int nd, size_t size, void** vv);
-void raid_par1_sse2(int nd, size_t size, void** vv);
-void raid_par2_int32(int nd, size_t size, void** vv);
-void raid_par2_int64(int nd, size_t size, void** vv);
-void raid_par2_sse2(int nd, size_t size, void** vv);
-void raid_par2_sse2ext(int nd, size_t size, void** vv);
-void raid_parz_int32(int nd, size_t size, void** vv);
-void raid_parz_int64(int nd, size_t size, void** vv);
-void raid_parz_sse2(int nd, size_t size, void** vv);
-void raid_parz_sse2ext(int nd, size_t size, void** vv);
-void raid_par3_int8(int nd, size_t size, void** vv);
-void raid_par3_ssse3(int nd, size_t size, void** vv);
-void raid_par3_ssse3ext(int nd, size_t size, void** vv);
-void raid_par4_int8(int nd, size_t size, void** vv);
-void raid_par4_ssse3(int nd, size_t size, void** vv);
-void raid_par4_ssse3ext(int nd, size_t size, void** vv);
-void raid_par5_int8(int nd, size_t size, void** vv);
-void raid_par5_ssse3(int nd, size_t size, void** vv);
-void raid_par5_ssse3ext(int nd, size_t size, void** vv);
-void raid_par6_int8(int nd, size_t size, void** vv);
-void raid_par6_ssse3(int nd, size_t size, void** vv);
-void raid_par6_ssse3ext(int nd, size_t size, void** vv);
-void raid_rec1_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
-void raid_rec2_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
-void raid_recX_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
-void raid_rec1_ssse3(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
-void raid_rec2_ssse3(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
-void raid_recX_ssse3(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
-
 #endif
diff --git a/tables.h b/raid/tables.c
similarity index 99%
rename from tables.h
rename to raid/tables.c
index 68011a55..cb5bda37 100644
--- a/tables.h
+++ b/raid/tables.c
@@ -15,7 +15,9 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-const uint8_t  __attribute__((aligned(256))) gfmul[256][256] =
+#include "internal.h"
+
+const uint8_t __attribute__((aligned(256))) raid_gfmul[256][256] =
 {
 	{
 		0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
@@ -8723,7 +8725,7 @@ const uint8_t  __attribute__((aligned(256))) gfmul[256][256] =
 	},
 };
 
-const uint8_t __attribute__((aligned(256))) gfexp[256] =
+const uint8_t __attribute__((aligned(256))) raid_gfexp[256] =
 {
 	0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80,
 	0x1d, 0x3a, 0x74, 0xe8, 0xcd, 0x87, 0x13, 0x26,
@@ -8759,7 +8761,7 @@ const uint8_t __attribute__((aligned(256))) gfexp[256] =
 	0x1b, 0x36, 0x6c, 0xd8, 0xad, 0x47, 0x8e, 0x01,
 };
 
-const uint8_t __attribute__((aligned(256))) gfinv[256] =
+const uint8_t __attribute__((aligned(256))) raid_gfinv[256] =
 {
 	/* note that the first element is not significative */
 	0x00, 0x01, 0x8e, 0xf4, 0x47, 0xa7, 0x7a, 0xba,
@@ -8804,7 +8806,7 @@ const uint8_t __attribute__((aligned(256))) gfinv[256] =
  * 01 02 04 08 10 20 40 80 1d 3a 74 e8 cd 87 13 26 4c 98 2d 5a b4 75 ea c9 8f 03 06 0c 18 30 60 c0 9d 27 4e 9c 25 4a 94 35 6a d4 b5 77 ee c1 9f 23 46 8c 05 0a 14 28 50 a0 5d ba 69 d2 b9 6f de a1 5f be 61 c2 99 2f 5e bc 65 ca 89 0f 1e 3c 78 f0 fd e7 d3 bb 6b d6 b1 7f fe e1 df a3 5b b6 71 e2 d9 af 43 86 11 22 44 88 0d 1a 34 68 d0 bd 67 ce 81 1f 3e 7c f8 ed c7 93 3b 76 ec c5 97 33 66 cc 85 17 2e 5c b8 6d da a9 4f 9e 21 42 84 15 2a 54 a8 4d 9a 29 52 a4 55 aa 49 92 39 72 e4 d5 b7 73 e6 d1 bf 63 c6 91 3f 7e fc e5 d7 b3 7b f6 f1 ff e3 db ab 4b 96 31 62 c4 95 37 6e dc a5 57 ae 41 82 19 32 64 c8 8d 07 0e 1c 38 70 e0 dd a7 53 a6 51 a2 59 b2 79 f2 f9 ef c3 9b 2b 56 ac 45 8a 09 12 24 48 90 3d 7a f4 f5 f7 f3 fb eb cb 8b 0b 16 2c 58 b0 7d fa e9 cf 83 1b 36 6c 
  * 01 8e 47 ad d8 6c 36 1b 83 cf e9 fa 7d b0 58 2c 16 0b 8b cb eb fb f3 f7 f5 f4 7a 3d 90 48 24 12 09 8a 45 ac 56 2b 9b c3 ef f9 f2 79 b2 59 a2 51 a6 53 a7 dd e0 70 38 1c 0e 07 8d c8 64 32 19 82 41 ae 57 a5 dc 6e 37 95 c4 62 31 96 4b ab db e3 ff f1 f6 7b b3 d7 e5 fc 7e 3f 91 c6 63 bf d1 e6 73 b7 d5 e4 72 39 92 49 aa 55 a4 52 29 9a 4d a8 54 2a 15 84 42 21 9e 4f a9 da 6d b8 5c 2e 17 85 cc 66 33 97 c5 ec 76 3b 93 c7 ed f8 7c 3e 1f 81 ce 67 bd d0 68 34 1a 0d 88 44 22 11 86 43 af d9 e2 71 b6 5b a3 df e1 fe 7f b1 d6 6b bb d3 e7 fd f0 78 3c 1e 0f 89 ca 65 bc 5e 2f 99 c2 61 be 5f a1 de 6f b9 d2 69 ba 5d a0 50 28 14 0a 05 8c 46 23 9f c1 ee 77 b5 d4 6a 35 94 4a 25 9c 4e 27 9d c0 60 30 18 0c 06 03 8f c9 ea 75 b4 5a 2d 98 4c 26 13 87 cd e8 74 3a 1d 80 40 20 
  */
-const uint8_t  __attribute__((aligned(256))) gfvandermonde[3][256] =
+const uint8_t __attribute__((aligned(256))) raid_gfvandermonde[3][256] =
 {
 	{
 		0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01,
@@ -8921,7 +8923,7 @@ const uint8_t  __attribute__((aligned(256))) gfvandermonde[3][256] =
  * 01 97 7f 9c 7c 18 bd a2 58 1a da 74 70 a3 e5 47 29 07 f5 80 23 e9 fa 46 54 a0 99 95 53 9b 0b c7 09 c0 78 89 92 e3 0d b0 2a 8c fb 17 3f 26 65 87 27 5c 66 61 79 4d 32 b3 8d 52 e2 82 3d f9 c5 02 bc 4c 73 48 62 af ba 41 d9 c4 2f b1 33 b8 15 7d cf 3a a9 5f 84 6d 34 1b 44 94 72 81 42 be cc 4b 0a 6f 5a 22 36 b5 3c 9d 13 7e 08 dd d6 5e 04 fc 5b ec ef f1 6e 1e 77 24 e6 c6 aa cb fd 51 67 06 6a 4a 88 db b2 c2 5d 43 40 f7 50 a8 f2 7a 71 a4 d2 bf 31 90 19 9a 8e f6 c3 a6 e7 60 12 ee 2d de 38 e8 b7 98 c1 28 f3 05 96 63 d1 b9 14 9f 1d 83 68 75 ed 16 03 ce e4 df e0 10 ae 69 55 91 2e 4e fe 21 1f 9e e1 d5 cd ca f0 8b 2b c9 8a 93 bb 57 20 86 1c a1 4f 3e 25 d4 6c a5 6b a7 37 ff 39 35 0c f8 ea 56 45 8f 2c 59 ab 85 eb 49 0f dc d8 76 b6 f4 0e 11 b4 d0 30 d3 3b ad d7 
  * 01 2b 3f cf 73 2c d6 ed cb 74 15 78 8a c1 17 c9 89 68 21 ab 76 3b 4b 5a 6e 0e b9 d3 b6 3e 36 86 bf a2 a7 30 14 eb c7 2d 96 67 20 b5 9a e0 a8 c6 80 04 8d fe 75 5e 23 ca 8f 48 99 0d df 8e b8 70 29 9c 44 69 3d a5 c2 90 d2 1c 9b 02 1d 98 93 ec 84 e8 64 4c 3a 8b 97 f3 e5 c0 7d 26 c8 08 a0 62 82 55 f7 33 f6 51 63 4d 77 da fd c3 38 6d ee 09 47 a3 05 de a6 f1 22 25 6a 0c 81 b2 6b 58 d5 b3 fc fb 28 7f 07 dc 7a 9e d0 37 b4 e1 1a 24 03 ae 94 ba 88 2f ea 2e 8c 5b bb 79 d1 11 ff a4 19 3c 2a 4e 52 e3 95 bd 31 5d 35 4a 41 c4 db 42 c5 0b 49 1b 7c e4 b0 9d 45 f0 a9 61 57 06 d4 40 91 56 13 fa 87 ac 27 54 dd 59 1f 71 39 43 6c f9 be 4f f4 1e 32 cd e9 7e 7b 66 5f ef e7 6f 0a 60 d7 b7 83 92 e2 af 72 f8 b1 50 10 ce 18 53 a1 cc ad 12 34 0f f5 aa 16 e6 f2 d8 85 9f bc 
  */
-const uint8_t  __attribute__((aligned(256))) gfcauchy[6][256] =
+const uint8_t __attribute__((aligned(256))) raid_gfcauchy[6][256] =
 {
 	{
 		0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01,
@@ -9136,7 +9138,7 @@ const uint8_t  __attribute__((aligned(256))) gfcauchy[6][256] =
  * Indexes are [DISK][PARITY - 2][LH].
  * Where DISK is from 0 to 250, PARITY from 2 to 5, LH from 0 to 1.
  */
-const uint8_t  __attribute__((aligned(256))) gfcauchypshufb[251][4][2][16] =
+const uint8_t __attribute__((aligned(256))) raid_gfcauchypshufb[251][4][2][16] =
 {
 	{
 		{
@@ -13666,7 +13668,7 @@ const uint8_t  __attribute__((aligned(256))) gfcauchypshufb[251][4][2][16] =
  * Indexes are [MULTIPLER][LH].
  * Where MULTIPLER is from 0 to 255, LH from 0 to 1.
  */
-const uint8_t  __attribute__((aligned(256))) gfmulpshufb[256][2][16] =
+const uint8_t __attribute__((aligned(256))) raid_gfmulpshufb[256][2][16] =
 {
 	{
 		{ 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00 },
diff --git a/raid.c b/raid/x86.c
similarity index 58%
rename from raid.c
rename to raid/x86.c
index 07f5f32b..29f1922d 100644
--- a/raid.c
+++ b/raid/x86.c
@@ -15,269 +15,8 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-/*
- * Copyright notes:
- *
- * The RAID5 and RAID6 support was originally derived from the
- * H. Peter Anvin paper "The mathematics of RAID-6" [1] and the
- * libraid6 library, also by H. Peter Anvin, released with license
- * "GPLv2 or any later version" inside the Linux Kernel 2.6.38.
- *
- * This support was later completely rewritten (many times), but
- * the H. Peter Anvin's Copyright may still apply.
- *
- * The others RAID levels and the recovering based on matrix
- * inversion is original work implemented from scratch.
- */
-
-/*
- * The RAID5 and RAID6 support is implemented using the Galois Field
- * GF(2^8) with the primitive polynomial x^8 + x^4 + x^3 + x^2 + 1 
- * (285 decimal).
- *
- * The parity P and Q of a set of N disk Di with 0<=i<N, is computed
- * using the equations:
- *
- * P = sum(Di)
- * Q = sum(2^i * Di) with 0<=i<N
- *
- * This approach is the same used by the Linux Kernel RAID, and by ZFS
- * RAIDZ2, better described in the H. Peter Anvin paper "The mathematics 
- * of RAID-6" [1]. 
- *
- * To support triple parity, it was first evaluated and then dropped, an
- * extension of the same approach, with additional parity coefficients set
- * as powers of 2^-1, with equations:
- *
- * P = sum(Di)
- * Q = sum(2^i * Di)
- * R = sum(2^-i * Di) with 0<=i<N
- *
- * This approach works well for triple parity and it's very efficient,
- * because we can implement very fast parallel multiplications and
- * divisions by 2 in GF(2^8).
- *
- * It's also similar at the approach used by ZFS RAIDZ3, with the 
- * difference that ZFS uses powers of 4 instead of 2^-1.
- *
- * Unfortunately it doesn't work beyond triple parity, because whatever
- * value we choose to generate the power coefficients to compute other
- * parities, the resulting equations are not solvable for some
- * combinations of missing disks.
- *
- * This is expected, because the Vandermonde matrix used to compute the
- * parity has no guarantee to have all submatrices not singular
- * [2, Chap 11, Problem 7] and this is a requirement to have
- * a MDS (Maximum Distance Separable) code [2, Chap 11, Theorem 8].
- *
- * To overcome this limitation, we use a Cauchy matrix [3][4] to compute 
- * the parity. A Cauchy matrix has the property to have all the square
- * submatrices not singular, resulting in always solvable equations,
- * for any combination of missing disks.
- *
- * The problem of this approach is that it requires the use of
- * generic multiplications, and not only by 2 or 2^-1, potentially 
- * affecting badly the performance.
- *
- * Hopefully there is a method to implement parallel multiplications
- * using SSSE3 instructions [1][5]. Method competitive with the
- * computation of triple parity using power coefficients.
- *
- * Another important property of the Cauchy matrix is that we can setup
- * the first two rows with coeffients equal at the RAID5 and RAID6 approach
- * decribed, resulting in a compatible extension, and requiring SSSE3
- * instructions only if triple parity or beyond is used.
- * 
- * The matrix is also adjusted, multipling each row by a constant factor
- * to make the first column of all 1, to optimize the computation for
- * the first disk.
- *
- * This results in the matrix A[row,col] defined as:
- * 
- * 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01...
- * 01 02 04 08 10 20 40 80 1d 3a 74 e8 cd 87 13 26 4c 98 2d 5a b4 75...
- * 01 f5 d2 c4 9a 71 f1 7f fc 87 c1 c6 19 2f 40 55 3d ba 53 04 9c 61...
- * 01 bb a6 d7 c7 07 ce 82 4a 2f a5 9b b6 60 f1 ad e7 f4 06 d2 df 2e...
- * 01 97 7f 9c 7c 18 bd a2 58 1a da 74 70 a3 e5 47 29 07 f5 80 23 e9...
- * 01 2b 3f cf 73 2c d6 ed cb 74 15 78 8a c1 17 c9 89 68 21 ab 76 3b...
- * (see tables.h for the matrix with all the columns)
- *
- * This matrix supports 6 level of parity, one for each row, for up to 251
- * data disks, one for each column, with all the 377,342,351,231 square
- * submatrices not singular, verified also with brute-force.
- * 
- * This matrix can be extended to support any number of parities, just 
- * adding additional rows, and removing one column for each new row.
- * (see mktables.c for more details in how the matrix is generated)
- * 
- * In details, parity is computed as:
- *
- * P = sum(Di)
- * Q = sum(2^i *  Di)
- * R = sum(A[2,i] * Di)
- * S = sum(A[3,i] * Di)
- * T = sum(A[4,i] * Di)
- * U = sum(A[5,i] * Di) with 0<=i<N
- *
- * To recover from a failure of six disks at indexes x,y,z,h,v,w,
- * with 0<=x<y<z<h<v<w<N, we compute the parity of the available N-6 
- * disks as:
- *
- * Pa = sum(Di)
- * Qa = sum(2^i * Di)
- * Ra = sum(A[2,i] * Di)
- * Sa = sum(A[3,i] * Di)
- * Ta = sum(A[4,i] * Di)
- * Ua = sum(A[5,i] * Di) with 0<=i<N,i!=x,i!=y,i!=z,i!=h,i!=v,i!=w.
- *
- * And if we define:
- *
- * Pd = Pa + P
- * Qd = Qa + Q
- * Rd = Ra + R
- * Sd = Sa + S
- * Td = Ta + T
- * Ud = Ua + U
- *
- * we can sum these two sets of equations, obtaining:
- *
- * Pd =          Dx +          Dy +          Dz +          Dh +          Dv +          Dw
- * Qd =    2^x * Dx +    2^y * Dy +    2^z * Dz +    2^h * Dh +    2^v * Dv +    2^w * Dw
- * Rd = A[2,x] * Dx + A[2,y] * Dy + A[2,z] * Dz + A[2,h] * Dh + A[2,v] * Dv + A[2,w] * Dw
- * Sd = A[3,x] * Dx + A[3,y] * Dy + A[3,z] * Dz + A[3,h] * Dh + A[3,v] * Dv + A[3,w] * Dw
- * Td = A[4,x] * Dx + A[4,y] * Dy + A[4,z] * Dz + A[4,h] * Dh + A[4,v] * Dv + A[4,w] * Dw
- * Ud = A[5,x] * Dx + A[5,y] * Dy + A[5,z] * Dz + A[5,h] * Dh + A[5,v] * Dv + A[5,w] * Dw
- *
- * A linear system always solvable because the coefficients matrix is
- * always not singular due the properties of the matrix A[].
- *
- * Resulting speed in x64, with 8 data disks, using a stripe of 256 KiB,
- * for a Core i7-3740QM CPU @ 2.7GHz is:
- * 
- *             int8   int32   int64    sse2   sse2e   ssse3  ssse3e
- *   par1             11927   22075   36004
- *   par2              3378    5874   18235   19164
- *   par3       844                                    8814    9419
- *   par4       665                                    6836    7415
- *   par5       537                                    5388    5686
- *   par6       449                                    4307    4789
- *
- * Values are in MiB/s of data processed, not counting generated parity.
- *
- * You can replicate these results in your machine using the "snapraid -T"
- * command.
- *
- * For comparison, the triple parity computation using the power
- * coeffients "1,2,2^-1" is only a little faster than the one based on
- * the Cauchy matrix if SSSE3 is present.
- *
- *             int8   int32   int64    sse2   sse2e   ssse3  ssse3e
- *   parz              2112    3118    9589   10304
- *
- *
- * In conclusion, the use of power coefficients, and specifically powers
- * of 1,2,2^-1, is the best option to implement triple parity in CPUs 
- * without SSSE3.
- * But if a modern CPU with SSSE3 (or similar) is available, the Cauchy
- * matrix is the best option because it provides a fast and general 
- * approach working for any number of parities.
- *
- * References:
- * [1] Anvin, "The mathematics of RAID-6", 2004
- * [2] MacWilliams, Sloane, "The Theory of Error-Correcting Codes", 1977
- * [3] Blomer, "An XOR-Based Erasure-Resilient Coding Scheme", 1995
- * [4] Roth, "Introduction to Coding Theory", 2006
- * [5] Plank, "Screaming Fast Galois Field Arithmetic Using Intel SIMD Instructions", 2013
- */
-
-#include "portable.h"
-
-#include "raid.h"
-#include "cpu.h"
-#include "tables.h"
-
-/****************************************************************************/
-/* parity generation */
-
-/**
- * Generic assert() with inverted meaning.
- */
-#define bug_on(a) assert(!(a))
-
-/**
- * Dereference as uint8_t
- */
-#define v_8(p) (*(uint8_t*)&(p))
-
-/**
- * Dereference as uint32_t
- */
-#define v_32(p) (*(uint32_t*)&(p))
-
-/**
- * Dereference as uint64_t
- */
-#define v_64(p) (*(uint64_t*)&(p))
-
-/*
- * PAR1 (RAID5 with xor) 32bit C implementation
- */
-void raid_par1_int32(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	int d, l;
-	size_t i;
-
-	uint32_t p0;
-	uint32_t p1;
-
-	l = nd - 1;
-	p = v[nd];
-
-	for(i=0;i<size;i+=8) {
-		p0 = v_32(v[l][i]);
-		p1 = v_32(v[l][i+4]);
-		/* accessing disks in backward order because the buffers */
-		/* are also in backward order */
-		for(d=l-1;d>=0;--d) {
-			p0 ^= v_32(v[d][i]);
-			p1 ^= v_32(v[d][i+4]);
-		}
-		v_32(p[i]) = p0;
-		v_32(p[i+4]) = p1;
-	}
-}
-
-/*
- * PAR1 (RAID5 with xor) 64bit C implementation
- */
-void raid_par1_int64(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	int d, l;
-	size_t i;
-
-	uint64_t p0;
-	uint64_t p1;
-
-	l = nd - 1;
-	p = v[nd];
-
-	for(i=0;i<size;i+=16) {
-		p0 = v_64(v[l][i]);
-		p1 = v_64(v[l][i+8]);
-		/* accessing disks in backward order because the buffers */
-		/* are also in backward order */
-		for(d=l-1;d>=0;--d) {
-			p0 ^= v_64(v[d][i]);
-			p1 ^= v_64(v[d][i+8]);
-		}
-		v_64(p[i]) = p0;
-		v_64(p[i+8]) = p1;
-	}
-}
+#include "internal.h"
+#include "gf.h"
 
 #if defined(__i386__) || defined(__x86_64__)
 /*
@@ -325,147 +64,13 @@ void raid_par1_sse2(int nd, size_t size, void** vv)
 }
 #endif
 
-/**
- * Multiply each byte of a uint32 by 2 in the GF(2^8).
- */
-static __always_inline uint32_t x2_32(uint32_t v)
-{
-	uint32_t mask = v & 0x80808080U;
-	mask = (mask << 1) - (mask >> 7);
-	v = (v << 1) & 0xfefefefeU;
-	v ^= mask & 0x1d1d1d1dU;
-	return v;
-}
-
-/**
- * Multiply each byte of a uint64 by 2 in the GF(2^8).
- */
-static __always_inline uint64_t x2_64(uint64_t v)
-{
-	uint64_t mask = v & 0x8080808080808080ULL;
-	mask = (mask << 1) - (mask >> 7);
-	v = (v << 1) & 0xfefefefefefefefeULL;
-	v ^= mask & 0x1d1d1d1d1d1d1d1dULL;
-	return v;
-}
-
-/**
- * Divide each byte of a uint32 by 2 in the GF(2^8).
- */
-static __always_inline uint32_t d2_32(uint32_t v)
-{
-	uint32_t mask = v & 0x01010101U;
-	mask = (mask << 8) - mask;
-	v = (v >> 1) & 0x7f7f7f7fU;
-	v ^= mask & 0x8e8e8e8eU;
-	return v;
-}
-
-/**
- * Divide each byte of a uint64 by 2 in the GF(2^8).
- */
-static __always_inline uint64_t d2_64(uint64_t v)
-{
-	uint64_t mask = v & 0x0101010101010101ULL;
-	mask = (mask << 8) - mask;
-	v = (v >> 1) & 0x7f7f7f7f7f7f7f7fULL;
-	v ^= mask & 0x8e8e8e8e8e8e8e8eULL;
-	return v;
-}
-
-/*
- * PAR2 (RAID6 with powers of 2) 32bit C implementation
- */
-void raid_par2_int32(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	int d, l;
-	size_t i;
-
-	uint32_t d0, q0, p0;
-	uint32_t d1, q1, p1;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-
-	for(i=0;i<size;i+=8) {
-		q0 = p0 = v_32(v[l][i]);
-		q1 = p1 = v_32(v[l][i+4]);
-		for(d=l-1;d>=0;--d) {
-			d0 = v_32(v[d][i]);
-			d1 = v_32(v[d][i+4]);
-
-			p0 ^= d0;
-			p1 ^= d1;
-
-			q0 = x2_32(q0);
-			q1 = x2_32(q1);
-
-			q0 ^= d0;
-			q1 ^= d1;
-		}
-		v_32(p[i]) = p0;
-		v_32(p[i+4]) = p1;
-		v_32(q[i]) = q0;
-		v_32(q[i+4]) = q1;
-	}
-}
-
-/*
- * PAR2 (RAID6 with powers of 2) 64bit C implementation
- */
-void raid_par2_int64(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	int d, l;
-	size_t i;
-
-	uint64_t d0, q0, p0;
-	uint64_t d1, q1, p1;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-
-	for(i=0;i<size;i+=16) {
-		q0 = p0 = v_64(v[l][i]);
-		q1 = p1 = v_64(v[l][i+8]);
-		for(d=l-1;d>=0;--d) {
-			d0 = v_64(v[d][i]);
-			d1 = v_64(v[d][i+8]);
-
-			p0 ^= d0;
-			p1 ^= d1;
-
-			q0 = x2_64(q0);
-			q1 = x2_64(q1);
-
-			q0 ^= d0;
-			q1 ^= d1;
-		}
-		v_64(p[i]) = p0;
-		v_64(p[i+8]) = p1;
-		v_64(q[i]) = q0;
-		v_64(q[i+8]) = q1;
-	}
-}
-
 #if defined(__i386__) || defined(__x86_64__)
 static const struct gfconst16 {
 	uint8_t poly[16];
 	uint8_t low4[16];
-	uint8_t half[16];
-	uint8_t low7[16];
 } gfconst16  __attribute__((aligned(64))) = {
 	{ 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d },
 	{ 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f, 0x0f },
-	{ 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e },
-	{ 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f }
 };
 
 /*
@@ -597,291 +202,6 @@ void raid_par2_sse2ext(int nd, size_t size, void** vv)
 }
 #endif
 
-/*
- * PARz (triple parity with powers of 2^-1) 32bit C implementation
- */
-void raid_parz_int32(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	uint8_t* r;
-	int d, l;
-	size_t i;
-
-	uint32_t d0, r0, q0, p0;
-	uint32_t d1, r1, q1, p1;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-	r = v[nd+2];
-
-	for(i=0;i<size;i+=8) {
-		r0 = q0 = p0 = v_32(v[l][i]);
-		r1 = q1 = p1 = v_32(v[l][i+4]);
-		for(d=l-1;d>=0;--d) {
-			d0 = v_32(v[d][i]);
-			d1 = v_32(v[d][i+4]);
-
-			p0 ^= d0;
-			p1 ^= d1;
-
-			q0 = x2_32(q0);
-			q1 = x2_32(q1);
-
-			q0 ^= d0;
-			q1 ^= d1;
-
-			r0 = d2_32(r0);
-			r1 = d2_32(r1);
-
-			r0 ^= d0;
-			r1 ^= d1;
-		}
-		v_32(p[i]) = p0;
-		v_32(p[i+4]) = p1;
-		v_32(q[i]) = q0;
-		v_32(q[i+4]) = q1;
-		v_32(r[i]) = r0;
-		v_32(r[i+4]) = r1;
-	}
-}
-
-/*
- * PARz (triple parity with powers of 2^-1) 64bit C implementation
- */
-void raid_parz_int64(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	uint8_t* r;
-	int d, l;
-	size_t i;
-
-	uint64_t d0, r0, q0, p0;
-	uint64_t d1, r1, q1, p1;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-	r = v[nd+2];
-
-	for(i=0;i<size;i+=16) {
-		r0 = q0 = p0 = v_64(v[l][i]);
-		r1 = q1 = p1 = v_64(v[l][i+8]);
-		for(d=l-1;d>=0;--d) {
-			d0 = v_64(v[d][i]);
-			d1 = v_64(v[d][i+8]);
-
-			p0 ^= d0;
-			p1 ^= d1;
-
-			q0 = x2_64(q0);
-			q1 = x2_64(q1);
-
-			q0 ^= d0;
-			q1 ^= d1;
-
-			r0 = d2_64(r0);
-			r1 = d2_64(r1);
-
-			r0 ^= d0;
-			r1 ^= d1;
-		}
-		v_64(p[i]) = p0;
-		v_64(p[i+8]) = p1;
-		v_64(q[i]) = q0;
-		v_64(q[i+8]) = q1;
-		v_64(r[i]) = r0;
-		v_64(r[i+8]) = r1;
-	}
-}
-
-#if defined(__i386__) || defined(__x86_64__)
-/*
- * PARz (triple parity with powers of 2^-1) SSE2 implementation
- */
-void raid_parz_sse2(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	uint8_t* r;
-	int d, l;
-	size_t i;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-	r = v[nd+2];
-
-	asm volatile("movdqa %0,%%xmm7" : : "m" (gfconst16.poly[0]));
-	asm volatile("movdqa %0,%%xmm3" : : "m" (gfconst16.half[0]));
-	asm volatile("movdqa %0,%%xmm6" : : "m" (gfconst16.low7[0]));
-
-	for(i=0;i<size;i+=16) {
-		asm volatile("movdqa %0,%%xmm0" : : "m" (v[l][i]));
-		asm volatile("movdqa %xmm0,%xmm1");
-		asm volatile("movdqa %xmm0,%xmm2");
-		for(d=l-1;d>=0;--d) {
-			asm volatile("pxor %xmm4,%xmm4");
-			asm volatile("pcmpgtb %xmm1,%xmm4");
-			asm volatile("paddb %xmm1,%xmm1");
-			asm volatile("pand %xmm7,%xmm4");
-			asm volatile("pxor %xmm4,%xmm1");
-
-			asm volatile("movdqa %xmm2,%xmm4");
-			asm volatile("pxor %xmm5,%xmm5");
-			asm volatile("psllw $7,%xmm4");
-			asm volatile("psrlw $1,%xmm2");
-			asm volatile("pcmpgtb %xmm4,%xmm5");
-			asm volatile("pand %xmm6,%xmm2");
-			asm volatile("pand %xmm3,%xmm5");
-			asm volatile("pxor %xmm5,%xmm2");
-
-			asm volatile("movdqa %0,%%xmm4" : : "m" (v[d][i]));
-			asm volatile("pxor %xmm4,%xmm0");
-			asm volatile("pxor %xmm4,%xmm1");
-			asm volatile("pxor %xmm4,%xmm2");
-		}
-		asm volatile("movntdq %%xmm0,%0" : "=m" (p[i]));
-		asm volatile("movntdq %%xmm1,%0" : "=m" (q[i]));
-		asm volatile("movntdq %%xmm2,%0" : "=m" (r[i]));
-	}
-
-	asm volatile("sfence" : : : "memory");
-}
-#endif
-
-#if defined(__x86_64__)
-/*
- * PARz (triple parity with powers of 2^-1) SSE2 implementation
- *
- * Note that it uses 16 registers, meaning that x64 is required.
- */
-void raid_parz_sse2ext(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	uint8_t* r;
-	int d, l;
-	size_t i;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-	r = v[nd+2];
-
-	asm volatile("movdqa %0,%%xmm7" : : "m" (gfconst16.poly[0]));
-	asm volatile("movdqa %0,%%xmm3" : : "m" (gfconst16.half[0]));
-	asm volatile("movdqa %0,%%xmm11" : : "m" (gfconst16.low7[0]));
-
-	for(i=0;i<size;i+=32) {
-		asm volatile("movdqa %0,%%xmm0" : : "m" (v[l][i]));
-		asm volatile("movdqa %0,%%xmm8" : : "m" (v[l][i+16]));
-		asm volatile("movdqa %xmm0,%xmm1");
-		asm volatile("movdqa %xmm8,%xmm9");
-		asm volatile("movdqa %xmm0,%xmm2");
-		asm volatile("movdqa %xmm8,%xmm10");
-		for(d=l-1;d>=0;--d) {
-			asm volatile("movdqa %xmm2,%xmm6");
-			asm volatile("movdqa %xmm10,%xmm14");
-			asm volatile("pxor %xmm4,%xmm4");
-			asm volatile("pxor %xmm12,%xmm12");
-			asm volatile("pxor %xmm5,%xmm5");
-			asm volatile("pxor %xmm13,%xmm13");
-			asm volatile("psllw $7,%xmm6");
-			asm volatile("psllw $7,%xmm14");
-			asm volatile("psrlw $1,%xmm2");
-			asm volatile("psrlw $1,%xmm10");
-			asm volatile("pcmpgtb %xmm1,%xmm4");
-			asm volatile("pcmpgtb %xmm9,%xmm12");
-			asm volatile("pcmpgtb %xmm6,%xmm5");
-			asm volatile("pcmpgtb %xmm14,%xmm13");
-			asm volatile("paddb %xmm1,%xmm1");
-			asm volatile("paddb %xmm9,%xmm9");
-			asm volatile("pand %xmm11,%xmm2");
-			asm volatile("pand %xmm11,%xmm10");
-			asm volatile("pand %xmm7,%xmm4");
-			asm volatile("pand %xmm7,%xmm12");
-			asm volatile("pand %xmm3,%xmm5");
-			asm volatile("pand %xmm3,%xmm13");
-			asm volatile("pxor %xmm4,%xmm1");
-			asm volatile("pxor %xmm12,%xmm9");
-			asm volatile("pxor %xmm5,%xmm2");
-			asm volatile("pxor %xmm13,%xmm10");
-
-			asm volatile("movdqa %0,%%xmm4" : : "m" (v[d][i]));
-			asm volatile("movdqa %0,%%xmm12" : : "m" (v[d][i+16]));
-			asm volatile("pxor %xmm4,%xmm0");
-			asm volatile("pxor %xmm4,%xmm1");
-			asm volatile("pxor %xmm4,%xmm2");
-			asm volatile("pxor %xmm12,%xmm8");
-			asm volatile("pxor %xmm12,%xmm9");
-			asm volatile("pxor %xmm12,%xmm10");
-		}
-		asm volatile("movntdq %%xmm0,%0" : "=m" (p[i]));
-		asm volatile("movntdq %%xmm8,%0" : "=m" (p[i+16]));
-		asm volatile("movntdq %%xmm1,%0" : "=m" (q[i]));
-		asm volatile("movntdq %%xmm9,%0" : "=m" (q[i+16]));
-		asm volatile("movntdq %%xmm2,%0" : "=m" (r[i]));
-		asm volatile("movntdq %%xmm10,%0" : "=m" (r[i+16]));
-	}
-
-	asm volatile("sfence" : : : "memory");
-}
-#endif
-
-/*
- * PAR3 (triple parity with Cauchy matrix) 8bit C implementation
- * 
- * Note that instead of a generic multiplicationt table, likely resulting
- * in multiple cache misses, a precomputed table could be used. 
- * But this is only a kind of reference function, and we are not really
- * interested in speed.
- */
-void raid_par3_int8(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	uint8_t* r;
-	int d, l;
-	size_t i;
-
-	uint8_t d0, r0, q0, p0;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-	r = v[nd+2];
-
-	for(i=0;i<size;i+=1) {
-		p0 = q0 = r0 = 0;
-		for(d=l;d>0;--d) {
-			d0 = v_8(v[d][i]);
-
-			p0 ^= d0;
-			q0 ^= gfmul[d0][gfcauchy[1][d]];
-			r0 ^= gfmul[d0][gfcauchy[2][d]];
-		}
-
-		/* first disk with all coefficients at 1 */
-		d0 = v_8(v[0][i]);
-
-		p0 ^= d0;
-		q0 ^= d0;
-		r0 ^= d0;
-
-		v_8(p[i]) = p0;
-		v_8(q[i]) = q0;
-		v_8(r[i]) = r0;
-	}
-}
-
 #if defined(__i386__) || defined(__x86_64__)
 /*
  * PAR3 (triple parity with Cauchy matrix) SSSE3 implementation
@@ -1116,58 +436,6 @@ void raid_par3_ssse3ext(int nd, size_t size, void** vv)
 }
 #endif
 
-/*
- * PAR4 (quad parity with Cauchy matrix) 8bit C implementation
- * 
- * Note that instead of a generic multiplicationt table, likely resulting
- * in multiple cache misses, a precomputed table could be used. 
- * But this is only a kind of reference function, and we are not really
- * interested in speed.
- */
-void raid_par4_int8(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	uint8_t* r;
-	uint8_t* s;
-	int d, l;
-	size_t i;
-
-	uint8_t d0, s0, r0, q0, p0;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-	r = v[nd+2];
-	s = v[nd+3];
-
-	for(i=0;i<size;i+=1) {
-		p0 = q0 = r0 = s0 = 0;
-		for(d=l;d>0;--d) {
-			d0 = v_8(v[d][i]);
-
-			p0 ^= d0;
-			q0 ^= gfmul[d0][gfcauchy[1][d]];
-			r0 ^= gfmul[d0][gfcauchy[2][d]];
-			s0 ^= gfmul[d0][gfcauchy[3][d]];
-		}
-
-		/* first disk with all coefficients at 1 */
-		d0 = v_8(v[0][i]);
-
-		p0 ^= d0;
-		q0 ^= d0;
-		r0 ^= d0;
-		s0 ^= d0;
-
-		v_8(p[i]) = p0;
-		v_8(q[i]) = q0;
-		v_8(r[i]) = r0;
-		v_8(s[i]) = s0;
-	}
-}
-
 #if defined(__i386__) || defined(__x86_64__)
 /*
  * PAR4 (quad parity with Cauchy matrix) SSSE3 implementation
@@ -1453,63 +721,6 @@ void raid_par4_ssse3ext(int nd, size_t size, void** vv)
 }
 #endif
 
-/*
- * PAR5 (penta parity with Cauchy matrix) 8bit C implementation
- * 
- * Note that instead of a generic multiplicationt table, likely resulting
- * in multiple cache misses, a precomputed table could be used. 
- * But this is only a kind of reference function, and we are not really
- * interested in speed.
- */
-void raid_par5_int8(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	uint8_t* r;
-	uint8_t* s;
-	uint8_t* t;
-	int d, l;
-	size_t i;
-
-	uint8_t d0, t0, s0, r0, q0, p0;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-	r = v[nd+2];
-	s = v[nd+3];
-	t = v[nd+4];
-
-	for(i=0;i<size;i+=1) {
-		p0 = q0 = r0 = s0 = t0 = 0;
-		for(d=l;d>0;--d) {
-			d0 = v_8(v[d][i]);
-
-			p0 ^= d0;
-			q0 ^= gfmul[d0][gfcauchy[1][d]];
-			r0 ^= gfmul[d0][gfcauchy[2][d]];
-			s0 ^= gfmul[d0][gfcauchy[3][d]];
-			t0 ^= gfmul[d0][gfcauchy[4][d]];
-		}
-		
-		/* first disk with all coefficients at 1 */
-		d0 = v_8(v[0][i]);
-
-		p0 ^= d0;
-		q0 ^= d0;
-		r0 ^= d0;
-		s0 ^= d0;
-		t0 ^= d0;
-
-		v_8(p[i]) = p0;
-		v_8(q[i]) = q0;
-		v_8(r[i]) = r0;
-		v_8(s[i]) = s0;
-		v_8(t[i]) = t0;
-	}
-}
-
 #if defined(__i386__) || defined(__x86_64__)
 /*
  * PAR5 (penta parity with Cauchy matrix) SSSE3 implementation
@@ -1780,68 +991,6 @@ void raid_par5_ssse3ext(int nd, size_t size, void** vv)
 }
 #endif
 
-/*
- * PAR6 (hexa parity with Cauchy matrix) 8bit C implementation
- * 
- * Note that instead of a generic multiplicationt table, likely resulting
- * in multiple cache misses, a precomputed table could be used. 
- * But this is only a kind of reference function, and we are not really
- * interested in speed.
- */
-void raid_par6_int8(int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* q;
-	uint8_t* r;
-	uint8_t* s;
-	uint8_t* t;
-	uint8_t* u;
-	int d, l;
-	size_t i;
-
-	uint8_t d0, u0, t0, s0, r0, q0, p0;
-
-	l = nd - 1;
-	p = v[nd];
-	q = v[nd+1];
-	r = v[nd+2];
-	s = v[nd+3];
-	t = v[nd+4];
-	u = v[nd+5];
-
-	for(i=0;i<size;i+=1) {
-		p0 = q0 = r0 = s0 = t0 = u0 = 0;
-		for(d=l;d>0;--d) {
-			d0 = v_8(v[d][i]);
-
-			p0 ^= d0;
-			q0 ^= gfmul[d0][gfcauchy[1][d]];
-			r0 ^= gfmul[d0][gfcauchy[2][d]];
-			s0 ^= gfmul[d0][gfcauchy[3][d]];
-			t0 ^= gfmul[d0][gfcauchy[4][d]];
-			u0 ^= gfmul[d0][gfcauchy[5][d]];
-		}
-		
-		/* first disk with all coefficients at 1 */
-		d0 = v_8(v[0][i]);
-
-		p0 ^= d0;
-		q0 ^= d0;
-		r0 ^= d0;
-		s0 ^= d0;
-		t0 ^= d0;
-		u0 ^= d0;
-
-		v_8(p[i]) = p0;
-		v_8(q[i]) = q0;
-		v_8(r[i]) = r0;
-		v_8(s[i]) = s0;
-		v_8(t[i]) = t0;
-		v_8(u[i]) = u0;
-	}
-}
-
 #if defined(__i386__) || defined(__x86_64__)
 /*
  * PAR6 (hexa parity with Cauchy matrix) SSSE3 implementation
@@ -2150,469 +1299,6 @@ void raid_par6_ssse3ext(int nd, size_t size, void** vv)
 }
 #endif
 
-/* internal forwarder */
-static void (*raid_par_ptr[RAID_PARITY_MAX])(int nd, size_t size, void** vv);
-static void (*raid_par3_ptr)(int nd, size_t size, void** vv);
-static void (*raid_parz_ptr)(int nd, size_t size, void** vv);
-
-void raid_par(int np, int nd, size_t size, void** v)
-{
-	bug_on(np < 1 || np > RAID_PARITY_MAX);
-
-	raid_par_ptr[np - 1](nd, size, v);
-}
-
-/****************************************************************************/
-/* mode */
-
-/**
- * Generator matrix currently used.
- */
-static const uint8_t (*raid_gen)[256];
-
-void raid_mode(int mode)
-{
-	if (mode == RAID_MODE_VANDERMONDE) {
-		raid_par_ptr[2] = raid_parz_ptr;
-		raid_gen = gfvandermonde;
-	} else {
-		raid_par_ptr[2] = raid_par3_ptr;
-		raid_gen = gfcauchy;
-	}
-}
-
-/****************************************************************************/
-/* zero */
-
-/**
- * Buffer filled with 0 used in recovering.
- */
-static void* raid_empty_zero_page;
-
-void raid_zero(void* zero)
-{
-	raid_empty_zero_page = zero;
-}
-
-/****************************************************************************/
-/* recovering */
-
-/**
- * GF a*b.
- */
-static __always_inline uint8_t mul(uint8_t a, uint8_t b)
-{
-	return gfmul[a][b];
-}
-
-/**
- * GF 1/a.
- * Not defined for a == 0.
- */
-static __always_inline uint8_t inv(uint8_t v)
-{
-	bug_on(v == 0); /* division by zero */
-
-	return gfinv[v];
-}
-
-/**
- * GF 2^a.
- */
-static __always_inline uint8_t pow2(int v)
-{
-	bug_on(v < 0 || v > 254); /* invalid exponent */
-
-	return gfexp[v];
-}
-
-/**
- * Gets the multiplication table for a specified value.
- */
-static __always_inline const uint8_t* table(uint8_t v)
-{
-	return gfmul[v];
-}
-
-/**
- * Computes the parity without the missing data blocks
- * and store it in the buffers of such data blocks.
- *
- * This is the parity expressed as Pa,Qa,Ra,Sa,Ta,Ua
- * in the equations.
- *
- * Note that all the other parities not in the ip[] vector
- * are destroyed.
- */
-static void raid_delta_gen(int nr, const int* id, const int* ip, int nd, size_t size, void** v)
-{
-	void* p[RAID_PARITY_MAX];
-	void* pa[RAID_PARITY_MAX];
-	int i;
-
-	for(i=0;i<nr;++i) {
-		/* keep a copy of the parity buffer */
-		p[i] = v[nd+ip[i]];
-
-		/* buffer for missing data blocks */
-		pa[i] = v[id[i]];
-
-		/* set at zero the missing data blocks */
-		v[id[i]] = raid_empty_zero_page;
-
-		/* compute the parity over the missing data blocks */
-		v[nd+ip[i]] = pa[i];
-	}
-
-	/* recompute the minimal parity required */
-	raid_par(ip[nr - 1] + 1, nd, size, v);
-
-	for(i=0;i<nr;++i) {
-		/* restore disk buffers as before */
-		v[id[i]] = pa[i];
-
-		/* restore parity buffers as before */
-		v[nd+ip[i]] = p[i];
-	}
-}
-
-/**
- * Recover failure of one data block for PAR1.
- *
- * Starting from the equation:
- *
- * Pd = Dx
- *
- * and solving we get:
- *
- * Dx = Pd
- */
-static void raid_rec1_par1(const int* id, int nd, size_t size, void** v)
-{
-	void* p;
-	void* pa;
-
-	/* for PAR1 we can directly compute the missing block */
-	/* and we don't need to use the zero buffer */
-	p = v[nd];
-	pa = v[id[0]];
-
-	/* use the parity as missing data block */
-	v[id[0]] = p;
-
-	/* compute the parity over the missing data block */
-	v[nd] = pa;
-
-	/* compute */
-	raid_par(1, nd, size, v);
-
-	/* restore as before */
-	v[id[0]] = pa;
-	v[nd] = p;
-}
-
-/**
- * Recover failure of two data blocks for PAR2.
- *
- * Starting from the equations:
- *
- * Pd = Dx + Dy
- * Qd = 2^id[0] * Dx + 2^id[1] * Dy
- *
- * and solving we get:
- *
- *               1                     2^(-id[0])
- * Dy = ------------------- * Pd + ------------------- * Qd
- *      2^(id[1]-id[0]) + 1        2^(id[1]-id[0]) + 1
- *
- * Dx = Dy + Pd
- *
- * with conditions:
- *
- * 2^id[0] != 0
- * 2^(id[1]-id[0]) + 1 != 0
- *
- * That are always satisfied for any 0<=id[0]<id[1]<255.
- */
-static void raid_rec2_par2(const int* id, const int* ip, int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	size_t i;
-	uint8_t* p;
-	uint8_t* pa;
-	uint8_t* q;
-	uint8_t* qa;
-	const uint8_t* T[2];
-
-	/* get multiplication tables */
-	T[0] = table( inv(pow2(id[1]-id[0]) ^ 1) );
-	T[1] = table( inv(pow2(id[0]) ^ pow2(id[1])) );
-
-	/* compute delta parity */
-	raid_delta_gen(2, id, ip, nd, size, vv);
-
-	p = v[nd];
-	q = v[nd+1];
-	pa = v[id[0]];
-	qa = v[id[1]];
-
-	for(i=0;i<size;++i) {
-		/* delta */
-		uint8_t Pd = p[i] ^ pa[i];
-		uint8_t Qd = q[i] ^ qa[i];
-
-		/* reconstruct */
-		uint8_t Dy = T[0][Pd] ^ T[1][Qd];
-		uint8_t Dx = Pd ^ Dy;
-
-		/* set */
-		pa[i] = Dx;
-		qa[i] = Dy;
-	}
-}
-
-/****************************************************************************/
-/* matrix recovering */
-
-/**
- * Gets the generator matrix coefficient for parity 'p' and disk 'd'.
- */
-static __always_inline uint8_t A(int p, int d)
-{
-	return raid_gen[p][d];
-}
-
-/**
- * Inverts the square matrix M of size nxn into V.
- * We use Gauss elimination to invert.
- */
-static void raid_invert(uint8_t* M, uint8_t* V, int n)
-{
-	int i,j,k;
-
-	/* set the identity matrix in V */
-	for(i=0;i<n;++i) {
-		for(j=0;j<n;++j) {
-			V[i*n+j] = i == j;
-		}
-	}
-
-	/* for each element in the diagonal */
-	for(k=0;k<n;++k) {
-		uint8_t f;
-
-		/* the diagonal element cannot be 0 because */
-		/* we are inverting matrices with all the square submatrices */
-		/* not singular */
-		bug_on(M[k*n+k] == 0);
-
-		/* make the diagonal element to be 1 */
-		f = inv(M[k*n+k]);
-		for(j=0;j<n;++j) {
-			M[k*n+j] = mul(f, M[k*n+j]);
-			V[k*n+j] = mul(f, V[k*n+j]);
-		}
-
-		/* make all the elements over and under the diagonal to be 0 */
-		for(i=0;i<n;++i) {
-			if (i == k)
-				continue;
-			f = M[i*n+k];
-			for(j=0;j<n;++j) {
-				M[i*n+j] ^= mul(f, M[k*n+j]);
-				V[i*n+j] ^= mul(f, V[k*n+j]);
-			}
-		}
-	}
-}
-
-/**
- * Recover failure of one data block at index id[0] using parity at index ip[0] for any RAID level.
- *
- * Starting from the equation:
- *
- * Pd = A[ip[0],id[0]] * Dx
- *
- * and solving we get:
- *
- * Dx = A[ip[0],id[0]]^-1 * Pd
- */
-void raid_rec1_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* pa;
-	const uint8_t* T;
-	uint8_t G;
-	uint8_t V;
-	size_t i;
-
-	(void)nr; /* unused, it's always 1 */
-
-	/* if it's RAID5 uses the dedicated and faster function */
-	if (ip[0] == 0) {
-		raid_rec1_par1(id, nd, size, vv);
-		return;
-	}
-
-	/* setup the coefficients matrix */
-	G = A(ip[0],id[0]);
-
-	/* invert it to solve the system of linear equations */
-	V = inv(G);
-
-	/* get multiplication tables */
-	T = table(V);
-
-	/* compute delta parity */
-	raid_delta_gen(1, id, ip, nd, size, vv);
-
-	p = v[nd+ip[0]];
-	pa = v[id[0]];
-
-	for(i=0;i<size;++i) {
-		/* delta */
-		uint8_t Pd = p[i] ^ pa[i];
-
-		/* reconstruct */
-		pa[i] = T[Pd];
-	}
-}
-
-/**
- * Recover failure of two data blocks at indexes id[0],id[1] using parity at indexes ip[0],ip[1] for any RAID level.
- *
- * Starting from the equations:
- *
- * Pd = A[ip[0],id[0]] * Dx + A[ip[0],id[1]] * Dy
- * Qd = A[ip[1],id[0]] * Dx + A[ip[1],id[1]] * Dy
- *
- * we solve inverting the coefficients matrix.
- */
-void raid_rec2_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p;
-	uint8_t* pa;
-	uint8_t* q;
-	uint8_t* qa;
-	const int N = 2;
-	const uint8_t* T[N][N];
-	uint8_t G[N*N];
-	uint8_t V[N*N];
-	size_t i;
-	int j, k;
-
-	(void)nr; /* unused, it's always 2 */
-
-	/* if it's RAID6 uses the dedicated and faster function */
-	if (ip[0] == 0 && ip[1] == 1) {
-		raid_rec2_par2(id, ip, nd, size, vv);
-		return;
-	}
-
-	/* setup the coefficients matrix */
-	for(j=0;j<N;++j) {
-		for(k=0;k<N;++k) {
-			G[j*N+k] = A(ip[j],id[k]);
-		}
-	}
-
-	/* invert it to solve the system of linear equations */
-	raid_invert(G, V, N);
-
-	/* get multiplication tables */
-	for(j=0;j<N;++j) {
-		for(k=0;k<N;++k) {
-			T[j][k] = table( V[j*N+k] );
-		}
-	}
-
-	/* compute delta parity */
-	raid_delta_gen(2, id, ip, nd, size, vv);
-
-	p = v[nd+ip[0]];
-	q = v[nd+ip[1]];
-	pa = v[id[0]];
-	qa = v[id[1]];
-
-	for(i=0;i<size;++i) {
-		/* delta */
-		uint8_t Pd = p[i] ^ pa[i];
-		uint8_t Qd = q[i] ^ qa[i];
-
-		/* reconstruct */
-		pa[i] = T[0][0][Pd] ^ T[0][1][Qd];
-		qa[i] = T[1][0][Pd] ^ T[1][1][Qd];
-	}
-}
-
-/**
- * Recover failure of N data blocks at indexes id[N] using parity at indexes ip[N] for any RAID level.
- *
- * Starting from the N equations, with 0<=i<N :
- *
- * PD[i] = sum(A[ip[i],id[j]] * D[i]) 0<=j<N
- *
- * we solve inverting the coefficients matrix.
- *
- * Note that referring at previous equations you have:
- * PD[0] = Pd, PD[1] = Qd, PD[2] = Rd, ...
- * D[0] = Dx, D[1] = Dy, D[2] = Dz, ...
- */
-void raid_recX_int8(int nr, const int* id, const int* ip, int nd, size_t size, void** vv)
-{
-	uint8_t** v = (uint8_t**)vv;
-	uint8_t* p[RAID_PARITY_MAX];
-	uint8_t* pa[RAID_PARITY_MAX];
-	const uint8_t* T[RAID_PARITY_MAX][RAID_PARITY_MAX];
-	uint8_t G[RAID_PARITY_MAX*RAID_PARITY_MAX];
-	uint8_t V[RAID_PARITY_MAX*RAID_PARITY_MAX];
-	size_t i;
-	int j, k;
-
-	/* setup the coefficients matrix */
-	for(j=0;j<nr;++j) {
-		for(k=0;k<nr;++k) {
-			G[j*nr+k] = A(ip[j],id[k]);
-		}
-	}
-
-	/* invert it to solve the system of linear equations */
-	raid_invert(G, V, nr);
-
-	/* get multiplication tables */
-	for(j=0;j<nr;++j) {
-		for(k=0;k<nr;++k) {
-			T[j][k] = table( V[j*nr+k] );
-		}
-	}
-
-	/* compute delta parity */
-	raid_delta_gen(nr, id, ip, nd, size, vv);
-
-	for(j=0;j<nr;++j) {
-		p[j] = v[nd+ip[j]];
-		pa[j] = v[id[j]];
-	}
-
-	for(i=0;i<size;++i) {
-		uint8_t PD[RAID_PARITY_MAX];
-
-		/* delta */
-		for(j=0;j<nr;++j)
-			PD[j] = p[j][i] ^ pa[j][i];
-
-		/* reconstruct */
-		for(j=0;j<nr;++j) {
-			uint8_t b = 0;
-			for(k=0;k<nr;++k)
-				b ^= T[j][k][PD[k]];
-			pa[j][i] = b;
-		}
-	}
-}
-
 #if defined(__i386__) || defined(__x86_64__)
 /*
  * RAID recovering for one disk SSSE3 implementation
@@ -2841,7 +1527,6 @@ void raid_recX_ssse3(int nr, const int* id, const int* ip, int nd, size_t size,
 				asm volatile("pshufb %xmm5,%xmm3");
 				asm volatile("pxor   %xmm2,%xmm0");
 				asm volatile("pxor   %xmm3,%xmm1");
-				
 			}
 
 			asm volatile("pxor %xmm1,%xmm0");
@@ -2853,218 +1538,3 @@ void raid_recX_ssse3(int nr, const int* id, const int* ip, int nd, size_t size,
 }
 #endif
 
-/* internal forwarder */
-static void (*raid_rec_ptr[RAID_PARITY_MAX])(int nr, const int* id, const int* ip, int nd, size_t size, void** vv);
-
-void raid_rec(int nr, const int* id, const int* ip, int nd, size_t size, void** v)
-{
-	bug_on(nr > nd);
-
-	/* if failed data is present */
-	if (nr != 0) {
-		raid_rec_ptr[nr - 1](nr, id, ip, nd, size, v);
-	}
-}
-
-void raid_recpar(int nrd, const int* id, int nrp, int* ip, int np, int nd, size_t size, void** v)
-{
-	bug_on(nrd > nd);
-	bug_on(nrd + nrp > np);
-
-	/* if failed data is present */
-	if (nrd != 0) {
-		int iu[RAID_PARITY_MAX];
-		int i, j, k;
-
-		/* setup the vector of parities to use */
-		for(i=0,j=0,k=0;i<(int)np;++i) {
-			if (ip[j] == i) {
-				++j;
-			} else {
-				iu[k] = i;
-				++k;
-			}
-		}
-
-		/* recover the data */
-		raid_rec_ptr[nrd - 1](nrd, id, iu, nd, size, v);
-	}
-
-	/* recompute all the parities up to the last bad one */
-	if (nrp != 0) {
-		raid_par(ip[nrp - 1] + 1, nd, size, v);
-	}
-}
-
-/****************************************************************************/
-/* init/done */
-
-void raid_init(void)
-{
-	raid_par3_ptr = raid_par3_int8;
-	raid_par_ptr[3] = raid_par4_int8;
-	raid_par_ptr[4] = raid_par5_int8;
-	raid_par_ptr[5] = raid_par6_int8;
-
-	if (sizeof(void*) == 4) {
-		raid_par_ptr[0] = raid_par1_int32;
-		raid_par_ptr[1] = raid_par2_int32;
-		raid_parz_ptr = raid_parz_int32;
-	} else {
-		raid_par_ptr[0] = raid_par1_int64;
-		raid_par_ptr[1] = raid_par2_int64;
-		raid_parz_ptr = raid_parz_int64;
-	}
-
-	raid_rec_ptr[0] = raid_rec1_int8;
-	raid_rec_ptr[1] = raid_rec2_int8;
-	raid_rec_ptr[2] = raid_recX_int8;
-	raid_rec_ptr[3] = raid_recX_int8;
-	raid_rec_ptr[4] = raid_recX_int8;
-	raid_rec_ptr[5] = raid_recX_int8;
-
-#if defined(__i386__) || defined(__x86_64__)
-	if (raid_cpu_has_sse2()) {
-		raid_par_ptr[0] = raid_par1_sse2;
-#if defined(__x86_64__)
-		if (raid_cpu_has_slowextendedreg()) {
-			raid_par_ptr[1] = raid_par2_sse2;
-			raid_parz_ptr = raid_parz_sse2;
-		} else {
-			raid_par_ptr[1] = raid_par2_sse2ext;
-			raid_parz_ptr = raid_parz_sse2ext;
-		}
-#else
-		raid_par_ptr[1] = raid_par2_sse2;
-		raid_parz_ptr = raid_parz_sse2;
-#endif
-	}
-
-	if (raid_cpu_has_ssse3()) {
-#if defined(__x86_64__)
-		raid_par3_ptr = raid_par3_ssse3ext;
-		raid_par_ptr[3] = raid_par4_ssse3ext;
-		raid_par_ptr[4] = raid_par5_ssse3ext;
-		raid_par_ptr[5] = raid_par6_ssse3ext;
-#else
-		raid_par3_ptr = raid_par3_ssse3;
-		raid_par_ptr[3] = raid_par4_ssse3;
-		raid_par_ptr[4] = raid_par5_ssse3;
-		raid_par_ptr[5] = raid_par6_ssse3;
-#endif
-		raid_rec_ptr[0] = raid_rec1_ssse3;
-		raid_rec_ptr[1] = raid_rec2_ssse3;
-		raid_rec_ptr[2] = raid_recX_ssse3;
-		raid_rec_ptr[3] = raid_recX_ssse3;
-		raid_rec_ptr[4] = raid_recX_ssse3;
-		raid_rec_ptr[5] = raid_recX_ssse3;
-	}
-#endif
-
-	/* set the default mode */
-	raid_mode(RAID_MODE_CAUCHY);
-}
-
-static struct raid_func {
-	const char* name;
-	void* p;
-} RAID_FUNC[] = {
-	{ "int8", raid_par3_int8 },
-	{ "int8", raid_par4_int8 },
-	{ "int8", raid_par5_int8 },
-	{ "int8", raid_par6_int8 },
-	{ "int32", raid_par1_int32 },
-	{ "int64", raid_par1_int64 },
-	{ "int32", raid_par2_int32 },
-	{ "int64", raid_par2_int64 },
-	{ "int32", raid_parz_int32 },
-	{ "int64", raid_parz_int64 },
-	{ "int8", raid_rec1_int8 },
-	{ "int8", raid_rec2_int8 },
-	{ "int8", raid_recX_int8 },
-
-#if defined(__i386__) || defined(__x86_64__)
-	{ "sse2", raid_par1_sse2 },
-	{ "sse2", raid_par2_sse2 },
-	{ "sse2", raid_parz_sse2 },
-	{ "ssse3", raid_par3_ssse3 },
-	{ "ssse3", raid_par4_ssse3 },
-	{ "ssse3", raid_par5_ssse3 },
-	{ "ssse3", raid_par6_ssse3 },
-	{ "ssse3", raid_rec1_ssse3 },
-	{ "ssse3", raid_rec2_ssse3 },
-	{ "ssse3", raid_recX_ssse3 },
-#endif
-
-#if defined(__x86_64__)
-	{ "sse2e", raid_par2_sse2ext },
-	{ "sse2e", raid_parz_sse2ext },
-	{ "ssse3e", raid_par3_ssse3ext },
-	{ "ssse3e", raid_par4_ssse3ext },
-	{ "ssse3e", raid_par5_ssse3ext },
-	{ "ssse3e", raid_par6_ssse3ext },
-#endif
-	{ 0, 0 }
-};
-
-static const char* raid_tag(void* func)
-{
-	struct raid_func* i = RAID_FUNC;
-	while (i->name != 0) {
-		if (i->p == func)
-			return i->name;
-		++i;
-	}
-	return "unknown";
-}
-
-const char* raid_par1_tag(void)
-{
-	return raid_tag(raid_par_ptr[0]);
-}
-
-const char* raid_par2_tag(void)
-{
-	return raid_tag(raid_par_ptr[1]);
-}
-
-const char* raid_parz_tag(void)
-{
-	return raid_tag(raid_parz_ptr);
-}
-
-const char* raid_par3_tag(void)
-{
-	return raid_tag(raid_par_ptr[2]);
-}
-
-const char* raid_par4_tag(void)
-{
-	return raid_tag(raid_par_ptr[3]);
-}
-
-const char* raid_par5_tag(void)
-{
-	return raid_tag(raid_par_ptr[4]);
-}
-
-const char* raid_par6_tag(void)
-{
-	return raid_tag(raid_par_ptr[5]);
-}
-
-const char* raid_rec1_tag(void)
-{
-	return raid_tag(raid_rec_ptr[0]);
-}
-
-const char* raid_rec2_tag(void)
-{
-	return raid_tag(raid_rec_ptr[1]);
-}
-
-const char* raid_recX_tag(void)
-{
-	return raid_tag(raid_rec_ptr[2]);
-}
-
diff --git a/raid/x86z.c b/raid/x86z.c
new file mode 100644
index 00000000..bb96eb60
--- /dev/null
+++ b/raid/x86z.c
@@ -0,0 +1,165 @@
+/*
+ * Copyright (C) 2013 Andrea Mazzoleni
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "internal.h"
+
+#if defined(__i386__) || defined(__x86_64__)
+static const struct gfzconst16 {
+	uint8_t poly[16];
+	uint8_t half[16];
+	uint8_t low7[16];
+} gfzconst16  __attribute__((aligned(64))) = {
+	{ 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d, 0x1d },
+	{ 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e, 0x8e },
+	{ 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f, 0x7f }
+};
+
+/*
+ * PARz (triple parity with powers of 2^-1) SSE2 implementation
+ */
+void raid_parz_sse2(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	uint8_t* r;
+	int d, l;
+	size_t i;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+	r = v[nd+2];
+
+	asm volatile("movdqa %0,%%xmm7" : : "m" (gfzconst16.poly[0]));
+	asm volatile("movdqa %0,%%xmm3" : : "m" (gfzconst16.half[0]));
+	asm volatile("movdqa %0,%%xmm6" : : "m" (gfzconst16.low7[0]));
+
+	for(i=0;i<size;i+=16) {
+		asm volatile("movdqa %0,%%xmm0" : : "m" (v[l][i]));
+		asm volatile("movdqa %xmm0,%xmm1");
+		asm volatile("movdqa %xmm0,%xmm2");
+		for(d=l-1;d>=0;--d) {
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pcmpgtb %xmm1,%xmm4");
+			asm volatile("paddb %xmm1,%xmm1");
+			asm volatile("pand %xmm7,%xmm4");
+			asm volatile("pxor %xmm4,%xmm1");
+
+			asm volatile("movdqa %xmm2,%xmm4");
+			asm volatile("pxor %xmm5,%xmm5");
+			asm volatile("psllw $7,%xmm4");
+			asm volatile("psrlw $1,%xmm2");
+			asm volatile("pcmpgtb %xmm4,%xmm5");
+			asm volatile("pand %xmm6,%xmm2");
+			asm volatile("pand %xmm3,%xmm5");
+			asm volatile("pxor %xmm5,%xmm2");
+
+			asm volatile("movdqa %0,%%xmm4" : : "m" (v[d][i]));
+			asm volatile("pxor %xmm4,%xmm0");
+			asm volatile("pxor %xmm4,%xmm1");
+			asm volatile("pxor %xmm4,%xmm2");
+		}
+		asm volatile("movntdq %%xmm0,%0" : "=m" (p[i]));
+		asm volatile("movntdq %%xmm1,%0" : "=m" (q[i]));
+		asm volatile("movntdq %%xmm2,%0" : "=m" (r[i]));
+	}
+
+	asm volatile("sfence" : : : "memory");
+}
+#endif
+
+#if defined(__x86_64__)
+/*
+ * PARz (triple parity with powers of 2^-1) SSE2 implementation
+ *
+ * Note that it uses 16 registers, meaning that x64 is required.
+ */
+void raid_parz_sse2ext(int nd, size_t size, void** vv)
+{
+	uint8_t** v = (uint8_t**)vv;
+	uint8_t* p;
+	uint8_t* q;
+	uint8_t* r;
+	int d, l;
+	size_t i;
+
+	l = nd - 1;
+	p = v[nd];
+	q = v[nd+1];
+	r = v[nd+2];
+
+	asm volatile("movdqa %0,%%xmm7" : : "m" (gfzconst16.poly[0]));
+	asm volatile("movdqa %0,%%xmm3" : : "m" (gfzconst16.half[0]));
+	asm volatile("movdqa %0,%%xmm11" : : "m" (gfzconst16.low7[0]));
+
+	for(i=0;i<size;i+=32) {
+		asm volatile("movdqa %0,%%xmm0" : : "m" (v[l][i]));
+		asm volatile("movdqa %0,%%xmm8" : : "m" (v[l][i+16]));
+		asm volatile("movdqa %xmm0,%xmm1");
+		asm volatile("movdqa %xmm8,%xmm9");
+		asm volatile("movdqa %xmm0,%xmm2");
+		asm volatile("movdqa %xmm8,%xmm10");
+		for(d=l-1;d>=0;--d) {
+			asm volatile("movdqa %xmm2,%xmm6");
+			asm volatile("movdqa %xmm10,%xmm14");
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pxor %xmm12,%xmm12");
+			asm volatile("pxor %xmm5,%xmm5");
+			asm volatile("pxor %xmm13,%xmm13");
+			asm volatile("psllw $7,%xmm6");
+			asm volatile("psllw $7,%xmm14");
+			asm volatile("psrlw $1,%xmm2");
+			asm volatile("psrlw $1,%xmm10");
+			asm volatile("pcmpgtb %xmm1,%xmm4");
+			asm volatile("pcmpgtb %xmm9,%xmm12");
+			asm volatile("pcmpgtb %xmm6,%xmm5");
+			asm volatile("pcmpgtb %xmm14,%xmm13");
+			asm volatile("paddb %xmm1,%xmm1");
+			asm volatile("paddb %xmm9,%xmm9");
+			asm volatile("pand %xmm11,%xmm2");
+			asm volatile("pand %xmm11,%xmm10");
+			asm volatile("pand %xmm7,%xmm4");
+			asm volatile("pand %xmm7,%xmm12");
+			asm volatile("pand %xmm3,%xmm5");
+			asm volatile("pand %xmm3,%xmm13");
+			asm volatile("pxor %xmm4,%xmm1");
+			asm volatile("pxor %xmm12,%xmm9");
+			asm volatile("pxor %xmm5,%xmm2");
+			asm volatile("pxor %xmm13,%xmm10");
+
+			asm volatile("movdqa %0,%%xmm4" : : "m" (v[d][i]));
+			asm volatile("movdqa %0,%%xmm12" : : "m" (v[d][i+16]));
+			asm volatile("pxor %xmm4,%xmm0");
+			asm volatile("pxor %xmm4,%xmm1");
+			asm volatile("pxor %xmm4,%xmm2");
+			asm volatile("pxor %xmm12,%xmm8");
+			asm volatile("pxor %xmm12,%xmm9");
+			asm volatile("pxor %xmm12,%xmm10");
+		}
+		asm volatile("movntdq %%xmm0,%0" : "=m" (p[i]));
+		asm volatile("movntdq %%xmm8,%0" : "=m" (p[i+16]));
+		asm volatile("movntdq %%xmm1,%0" : "=m" (q[i]));
+		asm volatile("movntdq %%xmm9,%0" : "=m" (q[i+16]));
+		asm volatile("movntdq %%xmm2,%0" : "=m" (r[i]));
+		asm volatile("movntdq %%xmm10,%0" : "=m" (r[i+16]));
+	}
+
+	asm volatile("sfence" : : : "memory");
+}
+#endif
+
diff --git a/rehash.c b/rehash.c
index 261d7e2b..199d19cc 100644
--- a/rehash.c
+++ b/rehash.c
@@ -23,7 +23,7 @@
 #include "state.h"
 #include "parity.h"
 #include "handle.h"
-#include "raid.h"
+#include "raid/raid.h"
 
 /****************************************************************************/
 /* rehash */
diff --git a/scrub.c b/scrub.c
index 646dae39..b067d9a7 100644
--- a/scrub.c
+++ b/scrub.c
@@ -22,7 +22,7 @@
 #include "state.h"
 #include "parity.h"
 #include "handle.h"
-#include "raid.h"
+#include "raid/raid.h"
 
 /****************************************************************************/
 /* scrub */
diff --git a/selftest.c b/selftest.c
index aea71089..d509fd60 100644
--- a/selftest.c
+++ b/selftest.c
@@ -19,8 +19,9 @@
 
 #include "snapraid.h"
 #include "util.h"
-#include "raid.h"
-#include "cpu.h"
+#include "raid/raid.h"
+#include "raid/cpu.h"
+#include "raid/internal.h"
 #include "elem.h"
 #include "state.h"
 #include "combo.h"
diff --git a/snapraid.c b/snapraid.c
index 8d3f3ca0..e3174086 100644
--- a/snapraid.c
+++ b/snapraid.c
@@ -22,7 +22,7 @@
 #include "elem.h"
 #include "import.h"
 #include "state.h"
-#include "raid.h"
+#include "raid/raid.h"
 
 /****************************************************************************/
 /* main */
diff --git a/speed.c b/speed.c
index 1faf73fd..839ae53e 100644
--- a/speed.c
+++ b/speed.c
@@ -19,8 +19,9 @@
 
 #include "snapraid.h"
 #include "util.h"
-#include "raid.h"
-#include "cpu.h"
+#include "raid/raid.h"
+#include "raid/cpu.h"
+#include "raid/internal.h"
 #include "state.h"
 
 /**
diff --git a/state.c b/state.c
index 89d327d7..dee75b03 100644
--- a/state.c
+++ b/state.c
@@ -22,8 +22,8 @@
 #include "state.h"
 #include "util.h"
 #include "parity.h"
-#include "cpu.h"
-#include "raid.h"
+#include "raid/raid.h"
+#include "raid/cpu.h"
 
 const char* lev_name(unsigned l)
 {
diff --git a/status.c b/status.c
index 6d65ddcc..06ec1f71 100644
--- a/status.c
+++ b/status.c
@@ -22,7 +22,7 @@
 #include "state.h"
 #include "parity.h"
 #include "handle.h"
-#include "raid.h"
+#include "raid/raid.h"
 
 /****************************************************************************/
 /* status */
diff --git a/sync.c b/sync.c
index 9634326f..32e32151 100644
--- a/sync.c
+++ b/sync.c
@@ -22,7 +22,7 @@
 #include "state.h"
 #include "parity.h"
 #include "handle.h"
-#include "raid.h"
+#include "raid/raid.h"
 
 /****************************************************************************/
 /* sync */
diff --git a/util.c b/util.c
index f0efa796..00e08ae3 100644
--- a/util.c
+++ b/util.c
@@ -18,7 +18,7 @@
 #include "portable.h"
 
 #include "util.h"
-#include "cpu.h"
+#include "raid/cpu.h"
 
 /****************************************************************************/
 /* hex conversion table */