From f8303d29630c4f1b23b9e6e315284fe851dd9811 Mon Sep 17 00:00:00 2001
From: Andrea Mazzoleni <amadvance@gmail.com>
Date: Tue, 17 Sep 2013 23:05:44 +0200
Subject: [PATCH 1/4] Adds preliminary support for QuadParity.

---
 mktables.c |  16 +
 raid.c     | 983 ++++++++++++++++++++++++++++++++++++++++++++++++++++-
 raid.h     |  24 ++
 selftest.c | 264 +++++++++++++-
 speed.c    |  48 ++-
 state.c    |   5 +-
 state.h    |   6 +-
 tables.h   |  36 ++
 8 files changed, 1363 insertions(+), 19 deletions(-)

diff --git a/mktables.c b/mktables.c
index 15089baf..ee5d7c0b 100644
--- a/mktables.c
+++ b/mktables.c
@@ -98,6 +98,22 @@ int main(void)
 	}
 	printf("};\n\n");
 
+	/* 8^a */
+	printf("const unsigned char __attribute__((aligned(256))) gfexp8[256] =\n");
+	printf("{\n");
+	v = 1;
+	for(i=0;i<256;++i) {
+		if (i % 8 == 0)
+			printf("\t");
+		printf("0x%02x,", v);
+		v = gfmul(v, 8);
+		if (i % 8 == 7)
+			printf("\n");
+		else
+			printf(" ");
+	}
+	printf("};\n\n");
+
 	/* 1/a */
 	printf("const unsigned char __attribute__((aligned(256))) gfinv[256] =\n");
 	printf("{\n");
diff --git a/raid.c b/raid.c
index 93df72dc..bd470a85 100644
--- a/raid.c
+++ b/raid.c
@@ -708,7 +708,6 @@ void raidTP_sse2r1(unsigned char** buffer, unsigned diskmax, unsigned size)
 }
 #endif
 
-
 #if defined(__x86_64__)
 /*
  * RAIDTP unrolled-by-2 SSE2 implementation
@@ -794,6 +793,411 @@ void raidTP_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 }
 #endif
 
+/*
+ * RAIDQP unrolled-by-2 C implementation
+ */
+void raidQP_int32r2(unsigned char** buffer, unsigned diskmax, unsigned size)
+{
+	unsigned char* p;
+	unsigned char* q;
+	unsigned char* r;
+	unsigned char* s;
+	int z, z0;
+	unsigned d;
+
+	uint32_t wd0, ws0, wr0, wq0, wp0, w10, w20;
+	uint32_t wd1, ws1, wr1, wq1, wp1, w11, w21;
+
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
+	r = buffer[diskmax+2];
+	s = buffer[diskmax+3];
+
+	for(d=0;d<size;d+=8) {
+		ws0 = wr0 = wq0 = wp0 = *(uint32_t *)&buffer[z0][d+0*4];
+		ws1 = wr1 = wq1 = wp1 = *(uint32_t *)&buffer[z0][d+1*4];
+		for (z=z0-1;z>= 0;--z) {
+			wd0 = *(uint32_t *)&buffer[z][d+0*4];
+			wd1 = *(uint32_t *)&buffer[z][d+1*4];
+
+			/* wp += wd */
+			wp0 ^= wd0;
+			wp1 ^= wd1;
+
+			/* w1 = wq * {02} */
+			w20 = MASK(wq0);
+			w21 = MASK(wq1);
+			w10 = SHLBYTE(wq0);
+			w11 = SHLBYTE(wq1);
+			w20 &= NBYTES(0x1d);
+			w21 &= NBYTES(0x1d);
+			w10 ^= w20;
+			w11 ^= w21;
+
+			/* wq = w1 + wd */
+			wq0 = w10 ^ wd0;
+			wq1 = w11 ^ wd1;
+
+			/* w1 = wr * {02} */
+			w20 = MASK(wr0);
+			w21 = MASK(wr1);
+			w10 = SHLBYTE(wr0);
+			w11 = SHLBYTE(wr1);
+			w20 &= NBYTES(0x1d);
+			w21 &= NBYTES(0x1d);
+			w10 ^= w20;
+			w11 ^= w21;
+
+			/* w1 = w1 * {02} */
+			w20 = MASK(w10);
+			w21 = MASK(w11);
+			w10 = SHLBYTE(w10);
+			w11 = SHLBYTE(w11);
+			w20 &= NBYTES(0x1d);
+			w21 &= NBYTES(0x1d);
+			w10 ^= w20;
+			w11 ^= w21;
+
+			/* wr = w1 + wd */
+			wr0 = w10 ^ wd0;
+			wr1 = w11 ^ wd1;
+
+			/* w1 = ws * {02} */
+			w20 = MASK(ws0);
+			w21 = MASK(ws1);
+			w10 = SHLBYTE(ws0);
+			w11 = SHLBYTE(ws1);
+			w20 &= NBYTES(0x1d);
+			w21 &= NBYTES(0x1d);
+			w10 ^= w20;
+			w11 ^= w21;
+
+			/* w1 = w1 * {02} */
+			w20 = MASK(w10);
+			w21 = MASK(w11);
+			w10 = SHLBYTE(w10);
+			w11 = SHLBYTE(w11);
+			w20 &= NBYTES(0x1d);
+			w21 &= NBYTES(0x1d);
+			w10 ^= w20;
+			w11 ^= w21;
+
+			/* w1 = w1 * {02} */
+			w20 = MASK(w10);
+			w21 = MASK(w11);
+			w10 = SHLBYTE(w10);
+			w11 = SHLBYTE(w11);
+			w20 &= NBYTES(0x1d);
+			w21 &= NBYTES(0x1d);
+			w10 ^= w20;
+			w11 ^= w21;
+
+			/* ws = w1 + wd */
+			ws0 = w10 ^ wd0;
+			ws1 = w11 ^ wd1;
+
+		}
+		*(uint32_t *)&p[d+4*0] = wp0;
+		*(uint32_t *)&p[d+4*1] = wp1;
+		*(uint32_t *)&q[d+4*0] = wq0;
+		*(uint32_t *)&q[d+4*1] = wq1;
+		*(uint32_t *)&r[d+4*0] = wr0;
+		*(uint32_t *)&r[d+4*1] = wr1;
+		*(uint32_t *)&s[d+4*0] = ws0;
+		*(uint32_t *)&s[d+4*1] = ws1;
+	}
+}
+
+#if defined(__i386__) || defined(__x86_64__)
+/*
+ * RAIDQP unrolled-by-1 MMX implementation
+ */
+void raidQP_mmxr1(unsigned char** buffer, unsigned diskmax, unsigned size)
+{
+	unsigned char* p;
+	unsigned char* q;
+	unsigned char* r;
+	unsigned char* s;
+	int z, z0;
+	unsigned d;
+
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
+	r = buffer[diskmax+2];
+	s = buffer[diskmax+3];
+
+	asm volatile("movq %0,%%mm7" : : "m" (raid6_mmx_constants.x1d));
+	asm volatile("pxor %mm4,%mm4");
+	asm volatile("pxor %mm5,%mm5");
+
+	for(d=0;d<size;d+=8) {
+		asm volatile("movq %0,%%mm0" : : "m" (buffer[z0][d])); /* P[0] */
+		asm volatile("movq %mm0,%mm1"); /* Q[0] */
+		asm volatile("movq %mm0,%mm2"); /* R[0] */
+		asm volatile("movq %mm0,%mm3"); /* S[0] */
+		for(z=z0-1;z>=0;--z) {
+			asm volatile("pcmpgtb %mm1,%mm4");
+			asm volatile("pcmpgtb %mm3,%mm5");
+			asm volatile("paddb %mm1,%mm1");
+			asm volatile("paddb %mm3,%mm3");
+			asm volatile("pand %mm7,%mm4");
+			asm volatile("pand %mm7,%mm5");
+			asm volatile("pxor %mm4,%mm1");
+			asm volatile("pxor %mm5,%mm3");
+
+			asm volatile("pxor %mm4,%mm4");
+			asm volatile("pxor %mm5,%mm5");
+
+			asm volatile("pcmpgtb %mm2,%mm4");
+			asm volatile("pcmpgtb %mm3,%mm5");
+			asm volatile("paddb %mm2,%mm2");
+			asm volatile("paddb %mm3,%mm3");
+			asm volatile("pand %mm7,%mm4");
+			asm volatile("pand %mm7,%mm5");
+			asm volatile("pxor %mm4,%mm2");
+			asm volatile("pxor %mm5,%mm3");
+
+			asm volatile("pxor %mm4,%mm4");
+			asm volatile("pxor %mm5,%mm5");
+
+			asm volatile("pcmpgtb %mm2,%mm4");
+			asm volatile("pcmpgtb %mm3,%mm5");
+			asm volatile("paddb %mm2,%mm2");
+			asm volatile("paddb %mm3,%mm3");
+			asm volatile("pand %mm7,%mm4");
+			asm volatile("pand %mm7,%mm5");
+			asm volatile("pxor %mm4,%mm2");
+			asm volatile("pxor %mm5,%mm3");
+
+			asm volatile("movq %0,%%mm6" : : "m" (buffer[z][d]));
+			asm volatile("pxor %mm6,%mm0");
+			asm volatile("pxor %mm6,%mm1");
+			asm volatile("pxor %mm6,%mm2");
+			asm volatile("pxor %mm6,%mm3");
+
+			asm volatile("pxor %mm4,%mm4");
+			asm volatile("pxor %mm5,%mm5");
+		}
+		asm volatile("movq %%mm0,%0" : "=m" (p[d]));
+		asm volatile("movq %%mm1,%0" : "=m" (q[d]));
+		asm volatile("movq %%mm2,%0" : "=m" (r[d]));
+		asm volatile("movq %%mm3,%0" : "=m" (s[d]));
+	}
+
+	asm volatile("emms" : : : "memory");
+}
+#endif
+
+#if defined(__i386__) || defined(__x86_64__)
+/*
+ * RAIDQP unrolled-by-1 SSE2 implementation
+ */
+void raidQP_sse2r1(unsigned char** buffer, unsigned diskmax, unsigned size)
+{
+	unsigned char* p;
+	unsigned char* q;
+	unsigned char* r;
+	unsigned char* s;
+	int z, z0;
+	unsigned d;
+
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
+	r = buffer[diskmax+2];
+	s = buffer[diskmax+3];
+
+	asm volatile("movdqa %0,%%xmm7" : : "m" (raid6_sse_constants.x1d[0]));
+	asm volatile("pxor %xmm4,%xmm4");
+	asm volatile("pxor %xmm5,%xmm5");
+
+	for(d=0;d<size;d+=16) {
+		asm volatile("movdqa %0,%%xmm0" : : "m" (buffer[z0][d])); /* P[0] */
+		asm volatile("movdqa %xmm0,%xmm1"); /* Q[0] */
+		asm volatile("movdqa %xmm0,%xmm2"); /* R[0] */
+		asm volatile("movdqa %xmm0,%xmm3"); /* S[0] */
+		for(z=z0-1;z>=0;--z) {
+			asm volatile("pcmpgtb %xmm1,%xmm4");
+			asm volatile("pcmpgtb %xmm3,%xmm5");
+			asm volatile("paddb %xmm1,%xmm1");
+			asm volatile("paddb %xmm3,%xmm3");
+			asm volatile("pand %xmm7,%xmm4");
+			asm volatile("pand %xmm7,%xmm5");
+			asm volatile("pxor %xmm4,%xmm1");
+			asm volatile("pxor %xmm5,%xmm3");
+
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pxor %xmm5,%xmm5");
+
+			asm volatile("pcmpgtb %xmm2,%xmm4");
+			asm volatile("pcmpgtb %xmm3,%xmm5");
+			asm volatile("paddb %xmm2,%xmm2");
+			asm volatile("paddb %xmm3,%xmm3");
+			asm volatile("pand %xmm7,%xmm4");
+			asm volatile("pand %xmm7,%xmm5");
+			asm volatile("pxor %xmm4,%xmm2");
+			asm volatile("pxor %xmm5,%xmm3");
+
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pxor %xmm5,%xmm5");
+
+			asm volatile("pcmpgtb %xmm2,%xmm4");
+			asm volatile("pcmpgtb %xmm3,%xmm5");
+			asm volatile("paddb %xmm2,%xmm2");
+			asm volatile("paddb %xmm3,%xmm3");
+			asm volatile("pand %xmm7,%xmm4");
+			asm volatile("pand %xmm7,%xmm5");
+			asm volatile("pxor %xmm4,%xmm2");
+			asm volatile("pxor %xmm5,%xmm3");
+
+			asm volatile("movdqa %0,%%xmm6" : : "m" (buffer[z][d]));
+			asm volatile("pxor %xmm6,%xmm0");
+			asm volatile("pxor %xmm6,%xmm1");
+			asm volatile("pxor %xmm6,%xmm2");
+			asm volatile("pxor %xmm6,%xmm3");
+
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pxor %xmm5,%xmm5");
+		}
+		asm volatile("movntdq %%xmm0,%0" : "=m" (p[d]));
+		asm volatile("movntdq %%xmm1,%0" : "=m" (q[d]));
+		asm volatile("movntdq %%xmm2,%0" : "=m" (r[d]));
+		asm volatile("movntdq %%xmm3,%0" : "=m" (s[d]));
+	}
+
+	asm volatile("sfence" : : : "memory");
+}
+#endif
+
+#if defined(__x86_64__)
+/*
+ * RAIDQP unrolled-by-2 SSE2 implementation
+ * Note that it uses all the 16 registers, meaning that x64 is required.
+ */
+void raidQP_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
+{
+	unsigned char* p;
+	unsigned char* q;
+	unsigned char* r;
+	unsigned char* s;
+	int z, z0;
+	unsigned d;
+
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
+	r = buffer[diskmax+2];
+	s = buffer[diskmax+3];
+
+	asm volatile("movdqa %0,%%xmm7" : : "m" (raid6_sse_constants.x1d[0]));
+	asm volatile("pxor %xmm4,%xmm4");
+	asm volatile("pxor %xmm5,%xmm5");
+	asm volatile("pxor %xmm12,%xmm12");
+	asm volatile("pxor %xmm13,%xmm13");
+
+	for(d=0;d<size;d+=32) {
+		asm volatile("movdqa %0,%%xmm0" : : "m" (buffer[z0][d])); /* P[0] */
+		asm volatile("movdqa %0,%%xmm8" : : "m" (buffer[z0][d+16])); /* P[1] */
+		asm volatile("movdqa %xmm0,%xmm1"); /* Q[0] */
+		asm volatile("movdqa %xmm8,%xmm9"); /* Q[1] */
+		asm volatile("movdqa %xmm0,%xmm2"); /* R[0] */
+		asm volatile("movdqa %xmm8,%xmm10"); /* R[1] */
+		asm volatile("movdqa %xmm0,%xmm3"); /* S[0] */
+		asm volatile("movdqa %xmm8,%xmm11"); /* S[1] */
+		for(z=z0-1;z>=0;--z) {
+			asm volatile("pcmpgtb %xmm1,%xmm4");
+			asm volatile("pcmpgtb %xmm3,%xmm5");
+			asm volatile("pcmpgtb %xmm9,%xmm12");
+			asm volatile("pcmpgtb %xmm11,%xmm13");
+			asm volatile("paddb %xmm1,%xmm1");
+			asm volatile("paddb %xmm3,%xmm3");
+			asm volatile("paddb %xmm9,%xmm9");
+			asm volatile("paddb %xmm11,%xmm11");
+			asm volatile("pand %xmm7,%xmm4");
+			asm volatile("pand %xmm7,%xmm5");
+			asm volatile("pand %xmm7,%xmm12");
+			asm volatile("pand %xmm7,%xmm13");
+			asm volatile("pxor %xmm4,%xmm1");
+			asm volatile("pxor %xmm5,%xmm3");
+			asm volatile("pxor %xmm12,%xmm9");
+			asm volatile("pxor %xmm13,%xmm11");
+
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pxor %xmm5,%xmm5");
+			asm volatile("pxor %xmm12,%xmm12");
+			asm volatile("pxor %xmm13,%xmm13");
+
+			asm volatile("pcmpgtb %xmm2,%xmm4");
+			asm volatile("pcmpgtb %xmm3,%xmm5");
+			asm volatile("pcmpgtb %xmm10,%xmm12");
+			asm volatile("pcmpgtb %xmm11,%xmm13");
+			asm volatile("paddb %xmm2,%xmm2");
+			asm volatile("paddb %xmm3,%xmm3");
+			asm volatile("paddb %xmm10,%xmm10");
+			asm volatile("paddb %xmm11,%xmm11");
+			asm volatile("pand %xmm7,%xmm4");
+			asm volatile("pand %xmm7,%xmm5");
+			asm volatile("pand %xmm7,%xmm12");
+			asm volatile("pand %xmm7,%xmm13");
+			asm volatile("pxor %xmm4,%xmm2");
+			asm volatile("pxor %xmm5,%xmm3");
+			asm volatile("pxor %xmm12,%xmm10");
+			asm volatile("pxor %xmm13,%xmm11");
+
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pxor %xmm5,%xmm5");
+			asm volatile("pxor %xmm12,%xmm12");
+			asm volatile("pxor %xmm13,%xmm13");
+
+			asm volatile("pcmpgtb %xmm2,%xmm4");
+			asm volatile("pcmpgtb %xmm3,%xmm5");
+			asm volatile("pcmpgtb %xmm10,%xmm12");
+			asm volatile("pcmpgtb %xmm11,%xmm13");
+			asm volatile("paddb %xmm2,%xmm2");
+			asm volatile("paddb %xmm3,%xmm3");
+			asm volatile("paddb %xmm10,%xmm10");
+			asm volatile("paddb %xmm11,%xmm11");
+			asm volatile("pand %xmm7,%xmm4");
+			asm volatile("pand %xmm7,%xmm5");
+			asm volatile("pand %xmm7,%xmm12");
+			asm volatile("pand %xmm7,%xmm13");
+			asm volatile("pxor %xmm4,%xmm2");
+			asm volatile("pxor %xmm5,%xmm3");
+			asm volatile("pxor %xmm12,%xmm10");
+			asm volatile("pxor %xmm13,%xmm11");
+
+			asm volatile("movdqa %0,%%xmm6" : : "m" (buffer[z][d]));
+			asm volatile("movdqa %0,%%xmm14" : : "m" (buffer[z][d+16]));
+			asm volatile("pxor %xmm6,%xmm0");
+			asm volatile("pxor %xmm6,%xmm1");
+			asm volatile("pxor %xmm6,%xmm2");
+			asm volatile("pxor %xmm6,%xmm3");
+			asm volatile("pxor %xmm14,%xmm8");
+			asm volatile("pxor %xmm14,%xmm9");
+			asm volatile("pxor %xmm14,%xmm10");
+			asm volatile("pxor %xmm14,%xmm11");
+
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pxor %xmm5,%xmm5");
+			asm volatile("pxor %xmm12,%xmm12");
+			asm volatile("pxor %xmm13,%xmm13");
+		}
+		asm volatile("movntdq %%xmm0,%0" : "=m" (p[d]));
+		asm volatile("movntdq %%xmm8,%0" : "=m" (p[d+16]));
+		asm volatile("movntdq %%xmm1,%0" : "=m" (q[d]));
+		asm volatile("movntdq %%xmm9,%0" : "=m" (q[d+16]));
+		asm volatile("movntdq %%xmm2,%0" : "=m" (r[d]));
+		asm volatile("movntdq %%xmm10,%0" : "=m" (r[d+16]));
+		asm volatile("movntdq %%xmm3,%0" : "=m" (s[d]));
+		asm volatile("movntdq %%xmm11,%0" : "=m" (s[d+16]));
+	}
+
+	asm volatile("sfence" : : : "memory");
+}
+#endif
+
 /****************************************************************************/
 /* generic computation */
 
@@ -801,6 +1205,7 @@ void raidTP_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 static void (*raid5_gen)(unsigned char** buffer, unsigned diskmax, unsigned size);
 static void (*raid6_gen)(unsigned char** buffer, unsigned diskmax, unsigned size);
 static void (*raidTP_gen)(unsigned char** buffer, unsigned diskmax, unsigned size);
+static void (*raidQP_gen)(unsigned char** buffer, unsigned diskmax, unsigned size);
 
 void raid_gen(unsigned level, unsigned char** buffer, unsigned diskmax, unsigned size)
 {
@@ -808,6 +1213,7 @@ void raid_gen(unsigned level, unsigned char** buffer, unsigned diskmax, unsigned
 	case 1 : raid5_gen(buffer, diskmax, size); break;
 	case 2 : raid6_gen(buffer, diskmax, size); break;
 	case 3 : raidTP_gen(buffer, diskmax, size); break;
+	case 4 : raidQP_gen(buffer, diskmax, size); break;
 	}
 }
 
@@ -914,45 +1320,68 @@ static inline unsigned char mul(unsigned char a, unsigned char b)
  */
 static inline unsigned char inv(unsigned char a)
 {
-#ifdef CHECKER
 	if (a == 0) {
 		fprintf(stderr, "GF division by zero\n");
 		exit(EXIT_FAILURE);
 	}
-#endif
+
 	return gfinv[a];
 }
 
 /**
  * GF 2^a.
- * Not defined for a == 255.
  */
 static inline unsigned char pow2(int v)
 {
-#ifdef CHECKER
 	if (v < 0 || v > 254) {
 		fprintf(stderr, "GF invalid exponent\n");
 		exit(EXIT_FAILURE);
 	}
-#endif
+
 	return gfexp2[v];
 }
-
 /**
  * GF 4^a.
- * Not defined for a == 255.
  */
 static inline unsigned char pow4(int v)
 {
-#ifdef CHECKER
 	if (v < 0 || v > 254) {
 		fprintf(stderr, "GF invalid exponent\n");
 		exit(EXIT_FAILURE);
 	}
-#endif
+
 	return gfexp4[v];
 }
 
+/**
+ * GF 8^a.
+ */
+static inline unsigned char pow8(int v)
+{
+	if (v < 0 || v > 254) {
+		fprintf(stderr, "GF invalid exponent\n");
+		exit(EXIT_FAILURE);
+	}
+
+	return gfexp8[v];
+}
+
+/**
+ * GF 2^b^a.
+ */
+static inline unsigned char pown(int b, int a)
+{
+	switch (b) {
+	case 0 : return 1;
+	case 1 : return pow2(a);
+	case 2 : return pow4(a);
+	case 3 : return pow8(a);
+	}
+
+	fprintf(stderr, "GF invalid base\n");
+	exit(EXIT_FAILURE);
+}
+
 /**
  * GF multiplication table.
  */
@@ -1480,6 +1909,536 @@ void raidTP_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size,
 	}
 }
 
+/****************************************************************************/
+/* matrix recovering  */
+
+/**
+ * Inverts the square matrix M of size nxn into V.
+ * We use Gauss elimination to invert.
+ */
+static void invert(unsigned char* M, unsigned char* V, unsigned n)
+{
+	unsigned 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) {
+		unsigned char f;
+
+		/* if the diagonal element is 0 */
+		if (M[k*n+k] == 0) {
+			/* find the next row with element at the same position not 0 */
+			for(i=k+1;i<n;++i)
+				if (M[i*n+k] != 0)
+					break;
+
+			/* if not found, the matrix is singular */
+			if (i == n) {
+				fprintf(stderr, "GF inversion of a singular matrix\n");
+				exit(EXIT_FAILURE);
+			}
+
+			/* sum the row at the present one to make the diagonal element different than 0 */
+			/* here usually the two rows are swapped, but summing is also ok */
+			for(j=0;j<n;++j) {
+				M[k*n+j] ^= M[i*n+j];
+				V[k*n+j] ^= V[i*n+j];
+			}
+		}
+
+		/* 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]);
+			}
+		}
+	}
+}
+
+/*
+ * Starting from the equation:
+ *
+ * Sd = 2^i^x * Dx
+ *
+ * and solving we get:
+ *
+ * Dx = 2^i^(-x) * Sd
+ *
+ * with conditions:
+ *
+ * 2^i^x != 0
+ *
+ * That are always satisfied for any 0<=x<y<255.
+ */
+void raid_recov_1data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int i, unsigned char* zero)
+{
+	unsigned char* p;
+	unsigned char* dp;
+	const unsigned char* pxmul; /* P multiplier to compute Dx */
+
+	/* if i is P use the fast one */
+	if (i == 0)
+		raid5_recov_data(dptrs, diskmax, size, x);
+
+	p = dptrs[diskmax+i];
+
+	/* compute syndrome with zero for the missing data page */
+	/* use the dead data page as temporary storage for delta p */
+	dp = dptrs[x];
+	dptrs[x] = zero;
+	dptrs[diskmax+i] = dp;
+
+	raid_gen(i + 1, dptrs, diskmax, size);
+
+	/* restore pointers */
+	dptrs[x] = dp;
+	dptrs[diskmax+i] = p;
+
+	/* select tables */
+	pxmul = table( inv(pown(i,x)) );
+
+	while (size--) {
+		/* delta */
+		unsigned char Pd = *p ^ *dp;
+
+		/* addends to reconstruct Dx */
+		unsigned char pxm = pxmul[Pd];
+
+		/* reconstruct Dx */
+		unsigned char Dx = pxm;
+
+		/* set */
+		*dp = Dx;
+
+		++p;
+		++dp;
+	}
+}
+
+/*
+ * Starting from the equations:
+ *
+ * Pd = 2^i^x * Dx + 2^i^y * Dy
+ * Qd = 2^j^x * Dx + 2^j^y * Dy
+ *
+ * we solve inverting the coefficients matrix.
+ */
+void raid_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, int i, int j, unsigned char* zero)
+{
+	unsigned char* p;
+	unsigned char* dp;
+	unsigned char* q;
+	unsigned char* dq;
+	const unsigned char* pxmul; /* P multiplier to compute Dx */
+	const unsigned char* qxmul; /* Q multiplier to compute Dx */
+	const unsigned char* pymul; /* P multiplier to compute Dy */
+	const unsigned char* qymul; /* Q multiplier to compute Dy */
+	unsigned char G[2*2];
+	unsigned char V[2*2];
+	unsigned n = 2;
+
+	p = dptrs[diskmax+i];
+	q = dptrs[diskmax+j];
+
+	/* compute syndrome with zero for the missing data page */
+	/* use the dead data page as temporary storage for delta p and q */
+	dp = dptrs[x];
+	dq = dptrs[y];
+	dptrs[x] = zero;
+	dptrs[y] = zero;
+	dptrs[diskmax+i] = dp;
+	dptrs[diskmax+j] = dq;
+
+	raid_gen(j + 1, dptrs, diskmax, size);
+
+	/* restore pointers */
+	dptrs[x] = dp;
+	dptrs[y] = dq;
+	dptrs[diskmax+i] = p;
+	dptrs[diskmax+j] = q;
+
+	/* setup the generator matrix */
+	G[0*n+0] = pown(i,x); /* row 1 for P */
+	G[0*n+1] = pown(i,y);
+	G[1*n+0] = pown(j,x); /* row 2 for Q */
+	G[1*n+1] = pown(j,y);
+
+	/* invert it to solve the system of linear equations */
+	invert(G, V, n);
+
+	/* select tables */
+	pxmul = table( V[0*n+0] );
+	qxmul = table( V[0*n+1] );
+	pymul = table( V[1*n+0] );
+	qymul = table( V[1*n+1] );
+
+	while (size--) {
+		/* delta */
+		unsigned char Pd = *p ^ *dp;
+		unsigned char Qd = *q ^ *dq;
+
+		/* addends to reconstruct Dy */
+		unsigned char pym = pymul[Pd];
+		unsigned char qym = qymul[Qd];
+
+		/* reconstruct Dy */
+		unsigned char Dy = pym ^ qym;
+
+		/* reconstruct Dx */
+		unsigned char Dx;
+
+		/* if i is P, take the fast way */
+		if (i == 0) {
+			Dx = Pd ^ Dy;
+		} else {
+			/* addends to reconstruct Dx */
+			unsigned char pxm = pxmul[Pd];
+			unsigned char qxm = qxmul[Qd];
+
+			/* reconstruct Dx */
+			Dx = pxm ^ qxm;
+		}
+
+		/* set */
+		*dp = Dx;
+		*dq = Dy;
+
+		++p;
+		++dp;
+		++q;
+		++dq;
+	}
+}
+
+/*
+ * Starting from the equations:
+ *
+ * Pd = 2^i^x * Dx + 2^i^y * Dy + 2^i^z * Dz
+ * Qd = 2^j^x * Dx + 2^j^y * Dy + 2^j^z * Dz
+ * Rd = 2^k^x * Dx + 2^k^y * Dy + 2^k^z * Dz
+ *
+ * we solve inverting the coefficients matrix.
+ */
+void raid_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, int z, int i, int j, int k, unsigned char* zero)
+{
+	unsigned char* p;
+	unsigned char* dp;
+	unsigned char* q;
+	unsigned char* dq;
+	unsigned char* r;
+	unsigned char* dr;
+	const unsigned char* pxmul; /* P multiplier to compute Dx */
+	const unsigned char* qxmul; /* Q multiplier to compute Dx */
+	const unsigned char* rxmul; /* R multiplier to compute Dx */
+	const unsigned char* pymul; /* P multiplier to compute Dy */
+	const unsigned char* qymul; /* Q multiplier to compute Dy */
+	const unsigned char* rymul; /* R multiplier to compute Dy */
+	const unsigned char* pzmul; /* P multiplier to compute Dz */
+	const unsigned char* qzmul; /* Q multiplier to compute Dz */
+	const unsigned char* rzmul; /* R multiplier to compute Dz */
+	unsigned char G[3*3];
+	unsigned char V[3*3];
+	unsigned n = 3;
+
+	p = dptrs[diskmax+i];
+	q = dptrs[diskmax+j];
+	r = dptrs[diskmax+k];
+
+	/* compute syndrome with zero for the missing data page */
+	/* use the dead data page as temporary storage for delta p, q and r */
+	dp = dptrs[x];
+	dq = dptrs[y];
+	dr = dptrs[z];
+	dptrs[x] = zero;
+	dptrs[y] = zero;
+	dptrs[z] = zero;
+	dptrs[diskmax+i] = dp;
+	dptrs[diskmax+j] = dq;
+	dptrs[diskmax+k] = dr;
+
+	raid_gen(k + 1, dptrs, diskmax, size);
+
+	/* restore pointers */
+	dptrs[x] = dp;
+	dptrs[y] = dq;
+	dptrs[z] = dr;
+	dptrs[diskmax+i] = p;
+	dptrs[diskmax+j] = q;
+	dptrs[diskmax+k] = r;
+
+	/* setup the generator matrix */
+	G[0*n+0] = pown(i,x); /* row 1 for P */
+	G[0*n+1] = pown(i,y);
+	G[0*n+2] = pown(i,z);
+	G[1*n+0] = pown(j,x); /* row 2 for Q */
+	G[1*n+1] = pown(j,y);
+	G[1*n+2] = pown(j,z);
+	G[2*n+0] = pown(k,x); /* row 3 for R */
+	G[2*n+1] = pown(k,y);
+	G[2*n+2] = pown(k,z);
+
+	/* invert it to solve the system of linear equations */
+	invert(G, V, n);
+
+	/* select tables */
+	pxmul = table( V[0*n+0] );
+	qxmul = table( V[0*n+1] );
+	rxmul = table( V[0*n+2] );
+	pymul = table( V[1*n+0] );
+	qymul = table( V[1*n+1] );
+	rymul = table( V[1*n+2] );
+	pzmul = table( V[2*n+0] );
+	qzmul = table( V[2*n+1] );
+	rzmul = table( V[2*n+2] );
+
+	while (size--) {
+		/* delta */
+		unsigned char Pd = *p ^ *dp;
+		unsigned char Qd = *q ^ *dq;
+		unsigned char Rd = *r ^ *dr;
+
+		/* addends to reconstruct Dz */
+		unsigned char pzm = pzmul[Pd];
+		unsigned char qzm = qzmul[Qd];
+		unsigned char rzm = rzmul[Rd];
+
+		/* reconstruct Dz */
+		unsigned char Dz = pzm ^ qzm ^ rzm;
+
+		/* addends to reconstruct Dy */
+		unsigned char pym = pymul[Pd];
+		unsigned char qym = qymul[Qd];
+		unsigned char rym = rymul[Rd];
+
+		/* reconstruct Dy */
+		unsigned char Dy = pym ^ qym ^ rym;
+
+		/* reconstruct Dx */
+		unsigned char Dx;
+
+		/* if i is P, take the fast way */
+		if (i == 0) {
+			Dx = Pd ^ Dy ^ Dz;
+		} else {
+			/* addends to reconstruct Dx */
+			unsigned char pxm = pxmul[Pd];
+			unsigned char qxm = qxmul[Qd];
+			unsigned char rxm = rxmul[Rd];
+
+			/* reconstruct Dx */
+			Dx = pxm ^ qxm ^ rxm;
+		}
+
+		/* set */
+		*dp = Dx;
+		*dq = Dy;
+		*dr = Dz;
+
+		++p;
+		++dp;
+		++q;
+		++dq;
+		++r;
+		++dr;
+	}
+}
+
+/*
+ * Starting from the equations:
+ *
+ * Pd = 2^i^x * Dx + 2^i^y * Dy + 2^i^z * Dz + 2^i^v * Dv
+ * Qd = 2^j^x * Dx + 2^j^y * Dy + 2^j^z * Dz + 2^j^v * Dv
+ * Rd = 2^k^x * Dx + 2^k^y * Dy + 2^k^z * Dz + 2^k^v * Dv
+ * Sd = 2^l^x * Dx + 2^l^y * Dy + 2^l^z * Dz + 2^l^v * Dv
+ *
+ * we solve inverting the coefficients matrix.
+ */
+void raid_recov_4data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, int z, int v, int i, int j, int k, int l, unsigned char* zero)
+{
+	unsigned char* p;
+	unsigned char* dp;
+	unsigned char* q;
+	unsigned char* dq;
+	unsigned char* r;
+	unsigned char* dr;
+	unsigned char* s;
+	unsigned char* ds;
+	const unsigned char* pxmul; /* P multiplier to compute Dx */
+	const unsigned char* qxmul; /* Q multiplier to compute Dx */
+	const unsigned char* rxmul; /* R multiplier to compute Dx */
+	const unsigned char* sxmul; /* S multiplier to compute Dx */
+	const unsigned char* pymul; /* P multiplier to compute Dy */
+	const unsigned char* qymul; /* Q multiplier to compute Dy */
+	const unsigned char* rymul; /* R multiplier to compute Dy */
+	const unsigned char* symul; /* S multiplier to compute Dy */
+	const unsigned char* pzmul; /* P multiplier to compute Dz */
+	const unsigned char* qzmul; /* Q multiplier to compute Dz */
+	const unsigned char* rzmul; /* R multiplier to compute Dz */
+	const unsigned char* szmul; /* S multiplier to compute Dz */
+	const unsigned char* pvmul; /* P multiplier to compute Dv */
+	const unsigned char* qvmul; /* Q multiplier to compute Dv */
+	const unsigned char* rvmul; /* R multiplier to compute Dv */
+	const unsigned char* svmul; /* S multiplier to compute Dv */
+	unsigned char G[4*4];
+	unsigned char V[4*4];
+	unsigned n = 4;
+
+	p = dptrs[diskmax+i];
+	q = dptrs[diskmax+j];
+	r = dptrs[diskmax+k];
+	s = dptrs[diskmax+l];
+
+	/* compute syndrome with zero for the missing data page */
+	/* use the dead data page as temporary storage for delta p, q, r and s */
+	dp = dptrs[x];
+	dq = dptrs[y];
+	dr = dptrs[z];
+	ds = dptrs[v];
+	dptrs[x] = zero;
+	dptrs[y] = zero;
+	dptrs[z] = zero;
+	dptrs[v] = zero;
+	dptrs[diskmax+i] = dp;
+	dptrs[diskmax+j] = dq;
+	dptrs[diskmax+k] = dr;
+	dptrs[diskmax+l] = ds;
+
+	raid_gen(l + 1, dptrs, diskmax, size);
+
+	/* restore pointers */
+	dptrs[x] = dp;
+	dptrs[y] = dq;
+	dptrs[z] = dr;
+	dptrs[v] = ds;
+	dptrs[diskmax+i] = p;
+	dptrs[diskmax+j] = q;
+	dptrs[diskmax+k] = r;
+	dptrs[diskmax+l] = s;
+
+	/* setup the generator matrix */
+	G[0*n+0] = pown(i,x); /* row 1 for P */
+	G[0*n+1] = pown(i,y);
+	G[0*n+2] = pown(i,z);
+	G[0*n+3] = pown(i,v);
+	G[1*n+0] = pown(j,x); /* row 2 for Q */
+	G[1*n+1] = pown(j,y);
+	G[1*n+2] = pown(j,z);
+	G[1*n+3] = pown(j,v);
+	G[2*n+0] = pown(k,x); /* row 3 for R */
+	G[2*n+1] = pown(k,y);
+	G[2*n+2] = pown(k,z);
+	G[2*n+3] = pown(k,v);
+	G[3*n+0] = pown(l,x); /* row 4 for S */
+	G[3*n+1] = pown(l,y);
+	G[3*n+2] = pown(l,z);
+	G[3*n+3] = pown(l,v);
+
+	/* invert it to solve the system of linear equations */
+	invert(G, V, n);
+
+	/* select tables */
+	pxmul = table( V[0*n+0] );
+	qxmul = table( V[0*n+1] );
+	rxmul = table( V[0*n+2] );
+	sxmul = table( V[0*n+3] );
+	pymul = table( V[1*n+0] );
+	qymul = table( V[1*n+1] );
+	rymul = table( V[1*n+2] );
+	symul = table( V[1*n+3] );
+	pzmul = table( V[2*n+0] );
+	qzmul = table( V[2*n+1] );
+	rzmul = table( V[2*n+2] );
+	szmul = table( V[2*n+3] );
+	pvmul = table( V[3*n+0] );
+	qvmul = table( V[3*n+1] );
+	rvmul = table( V[3*n+2] );
+	svmul = table( V[3*n+3] );
+
+	while (size--) {
+		/* delta */
+		unsigned char Pd = *p ^ *dp;
+		unsigned char Qd = *q ^ *dq;
+		unsigned char Rd = *r ^ *dr;
+		unsigned char Sd = *s ^ *ds;
+
+		/* addends to reconstruct Dv */
+		unsigned char pvm = pvmul[Pd];
+		unsigned char qvm = qvmul[Qd];
+		unsigned char rvm = rvmul[Rd];
+		unsigned char svm = svmul[Sd];
+
+		/* reconstruct Dv */
+		unsigned char Dv = pvm ^ qvm ^ rvm ^ svm;
+
+		/* addends to reconstruct Dz */
+		unsigned char pzm = pzmul[Pd];
+		unsigned char qzm = qzmul[Qd];
+		unsigned char rzm = rzmul[Rd];
+		unsigned char szm = szmul[Sd];
+
+		/* reconstruct Dz */
+		unsigned char Dz = pzm ^ qzm ^ rzm ^ szm;
+
+		/* addends to reconstruct Dy */
+		unsigned char pym = pymul[Pd];
+		unsigned char qym = qymul[Qd];
+		unsigned char rym = rymul[Rd];
+		unsigned char sym = symul[Sd];
+
+		/* reconstruct Dy */
+		unsigned char Dy = pym ^ qym ^ rym ^ sym;
+
+		/* reconstruct Dx */
+		unsigned char Dx;
+
+		/* if i is P, take the fast way */
+		if (i == 0) {
+			Dx = Pd ^ Dy ^ Dz ^ Dv;
+		} else {
+			/* addends to reconstruct Dx */
+			unsigned char pxm = pxmul[Pd];
+			unsigned char qxm = qxmul[Qd];
+			unsigned char rxm = rxmul[Rd];
+			unsigned char sxm = sxmul[Sd];
+
+			/* reconstruct Dx */
+			Dx = pxm ^ qxm ^ rxm ^ sxm;
+		}
+
+		/* set */
+		*dp = Dx;
+		*dq = Dy;
+		*dr = Dz;
+		*ds = Dv;
+
+		++p;
+		++dp;
+		++q;
+		++dq;
+		++r;
+		++dr;
+		++s;
+		++ds;
+	}
+}
+
 /****************************************************************************/
 /* init/done */
 
@@ -1488,21 +2447,25 @@ void raid_init(void)
 	raid5_gen = raid5_int32r2;
 	raid6_gen = raid6_int32r2;
 	raidTP_gen = raidTP_int32r2;
+	raidQP_gen = raidQP_int32r2;
 #if defined(__i386__) || defined(__x86_64__)
 	if (cpu_has_mmx()) {
 		raid5_gen = raid5_mmxr4;
 		raid6_gen = raid6_mmxr2;
 		raidTP_gen = raidTP_mmxr1;
+		raidQP_gen = raidQP_mmxr1;
 	}
 	if (cpu_has_sse2()) {
 #if defined(__x86_64__)
 		raid5_gen = raid5_sse2r8;
 		raid6_gen = raid6_sse2r4;
 		raidTP_gen = raidTP_sse2r2;
+		raidQP_gen = raidQP_sse2r2;
 #else
 		raid5_gen = raid5_sse2r4;
 		raid6_gen = raid6_sse2r2;
 		raidTP_gen = raidTP_sse2r1;
+		raidQP_gen = raidQP_sse2r1;
 #endif
 	}
 #endif
diff --git a/raid.h b/raid.h
index 77011a48..120efe68 100644
--- a/raid.h
+++ b/raid.h
@@ -65,6 +65,26 @@ void raidTP_recov_2dataq(unsigned char** dptrs, unsigned diskmax, unsigned size,
  */
 void raidTP_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, int failc, unsigned char* zero);
 
+/*
+ * Recover failure of one data block x using parity i.
+ */
+void raid_recov_1data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int i, unsigned char* zero);
+
+/*
+ * Recover failure of two data blocks x,y using parity i,j.
+ */
+void raid_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, int i, int j, unsigned char* zero);
+
+/*
+ * Recover failure of three data blocks x,y,z using parity i,j,k.
+ */
+void raid_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, int z, int i, int j, int k, unsigned char* zero);
+
+/*
+ * Recover failure of four data blocks x,y,z,v using parity i,j,k,l.
+ */
+void raid_recov_4data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, int z, int v, int i, int j, int k, int l, unsigned char* zero);
+
 /**
  * Initializes the RAID system.
  */
@@ -87,6 +107,10 @@ void raidTP_int32r2(unsigned char** dptr, unsigned diskmax, unsigned size);
 void raidTP_mmxr1(unsigned char** buffer, unsigned diskmax, unsigned size);
 void raidTP_sse2r1(unsigned char** dptr, unsigned diskmax, unsigned size);
 void raidTP_sse2r2(unsigned char** dptr, unsigned diskmax, unsigned size);
+void raidQP_int32r2(unsigned char** dptr, unsigned diskmax, unsigned size);
+void raidQP_mmxr1(unsigned char** buffer, unsigned diskmax, unsigned size);
+void raidQP_sse2r1(unsigned char** dptr, unsigned diskmax, unsigned size);
+void raidQP_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size);
 
 #endif
 
diff --git a/selftest.c b/selftest.c
index 412f0cb7..01dd4ab3 100644
--- a/selftest.c
+++ b/selftest.c
@@ -374,6 +374,263 @@ static void raidTPtest(unsigned diskmax, unsigned block_size)
 	free(buffer);
 }
 
+static void raidQPtest(unsigned diskmax, unsigned block_size)
+{
+	void* buffer_alloc;
+	unsigned char** buffer;
+	unsigned char* buffer_test0;
+	unsigned char* buffer_test1;
+	unsigned char* buffer_test2;
+	unsigned char* buffer_test3;
+	unsigned char* buffer_zero;
+	unsigned char* buffer_bad;
+	unsigned buffermax;
+	unsigned i, j, k, l;
+	unsigned n;
+
+	struct combo1 {
+		int ignore[3];
+		int use;
+	} COMBO1[4] = {
+		{ { 1, 2, 3 }, 0 },
+		{ { 0, 2, 3 }, 1 },
+		{ { 0, 1, 3 }, 2 },
+		{ { 0, 1, 2 }, 3 }
+	};
+
+	struct combo2 {
+		int ignore[2];
+		int use[2];
+	} COMBO2[6] = {
+		{ { 0, 1 }, { 2, 3 } },
+		{ { 0, 2 }, { 1, 3 } },
+		{ { 0, 3 }, { 1, 2 } },
+		{ { 1, 2 }, { 0, 3 } },
+		{ { 1, 3 }, { 0, 2 } },
+		{ { 2, 3 }, { 0, 1 } }
+	};
+
+	struct combo3 {
+		int ignore;
+		int use[3];
+	} COMBO3[4] = {
+		{ 0, { 1, 2, 3 } },
+		{ 1, { 0, 2, 3 } },
+		{ 2, { 0, 1, 3 } },
+		{ 3, { 0, 1, 2 } }
+	};
+
+	buffermax = diskmax + 4 + 4 + 1 + 1;
+
+	buffer = malloc_nofail_vector_align(buffermax, block_size, &buffer_alloc);
+
+	buffer_test0 = buffer[buffermax-6];
+	buffer_test1 = buffer[buffermax-5];
+	buffer_test2 = buffer[buffermax-4];
+	buffer_test3 = buffer[buffermax-3];
+	buffer_zero = buffer[buffermax-2];
+	buffer_bad = buffer[buffermax-1];
+	memset(buffer_zero, 0, block_size);
+	memset(buffer_bad, 0xA5, block_size);
+
+	/* fill with random */
+	for(i=0;i<diskmax;++i) {
+		for(j=0;j<block_size;++j)
+			buffer[i][j] = rand();
+	}
+
+	/* compute the parity */
+	raid_gen(4, buffer, diskmax, block_size);
+
+	/* check 4data */
+	for(i=0;i<diskmax;++i) {
+		for(j=i+1;j<diskmax;++j) {
+			for(k=j+1;k<diskmax;++k) {
+				for(l=k+1;l<diskmax;++l) {
+					unsigned char* a;
+					unsigned char* b;
+					unsigned char* c;
+					unsigned char* d;
+
+					/* save */
+					a = buffer[i];
+					b = buffer[j];
+					c = buffer[k];
+					d = buffer[l];
+
+					/* damage */
+					buffer[i] = buffer_test0;
+					buffer[j] = buffer_test1;
+					buffer[k] = buffer_test2;
+					buffer[l] = buffer_test3;
+
+					/* recover */
+					raid_recov_4data(buffer, diskmax, block_size, i, j, k, l, 0, 1, 2, 3, buffer_zero);
+
+					/* check */
+					if (memcmp(buffer_test0, a, block_size) != 0
+						|| memcmp(buffer_test1, b, block_size) != 0
+						|| memcmp(buffer_test2, c, block_size) != 0
+						|| memcmp(buffer_test3, d, block_size) != 0
+					) {
+						fprintf(stderr, "Failed RAIDQP test\n");
+						exit(EXIT_FAILURE);
+					}
+
+					/* restore */
+					buffer[i] = a;
+					buffer[j] = b;
+					buffer[k] = c;
+					buffer[l] = d;
+				}
+			}
+		}
+	}
+
+	/* check 3data */
+	for(i=0;i<diskmax;++i) {
+		for(j=i+1;j<diskmax;++j) {
+			for(k=j+1;k<diskmax;++k) {
+				for(n=0;n<4;++n) {
+					unsigned char* a;
+					unsigned char* b;
+					unsigned char* c;
+					unsigned char* p;
+
+					/* save */
+					a = buffer[i];
+					b = buffer[j];
+					c = buffer[k];
+
+					/* damage */
+					buffer[i] = buffer_test0;
+					buffer[j] = buffer_test1;
+					buffer[k] = buffer_test2;
+
+					/* ignore parity */
+					l = diskmax + COMBO3[n].ignore;
+					p = buffer[l];
+					buffer[l] = buffer_bad;
+
+					/* recover ignoring p */
+					raid_recov_3data(buffer, diskmax, block_size, i, j, k, COMBO3[n].use[0], COMBO3[n].use[1], COMBO3[n].use[2], buffer_zero);
+
+					/* check */
+					if (memcmp(buffer_test0, a, block_size) != 0
+						|| memcmp(buffer_test1, b, block_size) != 0
+						|| memcmp(buffer_test2, c, block_size) != 0
+					) {
+						fprintf(stderr, "Failed RAIDQP test 1\n");
+						exit(EXIT_FAILURE);
+					}
+
+					/* restore parity */
+					buffer[l] = p;
+
+					/* restore */
+					buffer[i] = a;
+					buffer[j] = b;
+					buffer[k] = c;
+				}
+			}
+		}
+	}
+
+	/* check 2data */
+	for(i=0;i<diskmax;++i) {
+		for(j=i+1;j<diskmax;++j) {
+			for(n=0;n<6;++n) {
+				unsigned char* a;
+				unsigned char* b;
+				unsigned char* p;
+				unsigned char* q;
+
+				/* save */
+				a = buffer[i];
+				b = buffer[j];
+
+				/* damage */
+				buffer[i] = buffer_test0;
+				buffer[j] = buffer_test1;
+
+				/* ignore parity */
+				k = diskmax + COMBO2[n].ignore[0];
+				l = diskmax + COMBO2[n].ignore[1];
+				p = buffer[k];
+				q = buffer[l];
+				buffer[k] = buffer_bad;
+				buffer[l] = buffer_bad;
+
+				/* recover */
+				raid_recov_2data(buffer, diskmax, block_size, i, j, COMBO2[n].use[0], COMBO2[n].use[1], buffer_zero);
+
+				/* check */
+				if (memcmp(buffer_test0, a, block_size) != 0
+					|| memcmp(buffer_test1, b, block_size) != 0
+				) {
+					fprintf(stderr, "Failed RAIDQP test 1\n");
+					exit(EXIT_FAILURE);
+				}
+
+				/* restore parity */
+				buffer[k] = p;
+				buffer[l] = q;
+
+				/* restore */
+				buffer[i] = a;
+				buffer[j] = b;
+			}
+		}
+	}
+
+	/* check 1data */
+	for(n=0;n<4;++n) {
+		for(i=0;i<diskmax;++i) {
+			unsigned char* a;
+			unsigned char* p;
+			unsigned char* q;
+			unsigned char* r;
+
+			/* save */
+			a = buffer[i];
+
+			/* damage */
+			buffer[i] = buffer_test0;
+
+			/* ignore parity */
+			j = diskmax + COMBO1[n].ignore[0];
+			k = diskmax + COMBO1[n].ignore[1];
+			l = diskmax + COMBO1[n].ignore[2];
+			p = buffer[j];
+			q = buffer[k];
+			r = buffer[l];
+			buffer[j] = buffer_bad;
+			buffer[k] = buffer_bad;
+			buffer[l] = buffer_bad;
+
+			/* recover */
+			raid_recov_1data(buffer, diskmax, block_size, i, COMBO1[n].use, buffer_zero);
+
+			/* check */
+			if (memcmp(buffer_test0, a, block_size) != 0) {
+				fprintf(stderr, "Failed RAIDQP test\n");
+				exit(EXIT_FAILURE);
+			}
+
+			/* restore parity */
+			buffer[j] = p;
+			buffer[k] = q;
+			buffer[l] = r;
+
+			/* restore */
+			buffer[i] = a;
+		}
+	}
+
+	free(buffer_alloc);
+	free(buffer);
+}
+
 #define HASH_TEST_MAX 512 /* tests are never longer than 512 bytes */
 
 static void hashtest(void)
@@ -489,8 +746,9 @@ void selftest(int gui)
 
 	hashtest();
 	crc32ctest();
-	raid5test(33, 2048);
-	raid6test(33, 2048);
-	raidTPtest(33, 2048);
+	raid5test(33, 1024);
+	raid6test(33, 1024);
+	raidTPtest(33, 1024);
+	raidQPtest(21, 1024);
 }
 
diff --git a/speed.c b/speed.c
index c634e3a1..1a39cfee 100644
--- a/speed.c
+++ b/speed.c
@@ -66,10 +66,10 @@ void speed(void)
 	void* buffer_alloc;
 	unsigned char** buffer;
 
-	buffer = malloc_nofail_vector_align(diskmax + 3, block_size, &buffer_alloc);
+	buffer = malloc_nofail_vector_align(diskmax + 4, block_size, &buffer_alloc);
 
 	/* initialize with fixed data */
-	for(i=0;i<diskmax+3;++i) {
+	for(i=0;i<diskmax+4;++i) {
 		memset(buffer[i], i, block_size);
 	}
 	for(i=0;i<HASH_SIZE;++i)
@@ -161,6 +161,18 @@ void speed(void)
 		printf("sse2x2");
 #else
 		printf("sse2x1");
+#endif
+	else if (cpu_has_mmx())
+		printf("mmxx1");
+	else
+		printf("int32x2");
+	printf("\n");
+	printf("Expected fastest RAIDQP is ");
+	if (cpu_has_sse2())
+#if defined(__x86_64__)
+		printf("sse2x2");
+#else
+		printf("sse2x1");
 #endif
 	else if (cpu_has_mmx())
 		printf("mmxx1");
@@ -325,6 +337,38 @@ void speed(void)
 	}
 #endif
 
+	SPEED_START {
+		raidQP_int32r2(buffer, diskmax, block_size);
+	} SPEED_STOP
+
+	printf("RAIDQP int32x2 %"PRIu64" [MB/s]\n", ds / dt);
+
+#if defined(__i386__) || defined(__x86_64__)
+	if (cpu_has_mmx()) {
+		SPEED_START {
+			raidQP_mmxr1(buffer, diskmax, block_size);
+		} SPEED_STOP
+
+		printf("RAIDQP mmxr1 %"PRIu64" [MB/s]\n", ds / dt);
+	}
+
+	if (cpu_has_sse2()) {
+		SPEED_START {
+			raidQP_sse2r1(buffer, diskmax, block_size);
+		} SPEED_STOP
+
+		printf("RAIDQP sse2x1 %"PRIu64" [MB/s]\n", ds / dt);
+
+#if defined(__x86_64__)
+		SPEED_START {
+			raidQP_sse2r2(buffer, diskmax, block_size);
+		} SPEED_STOP
+
+		printf("RAIDQP sse2x2 %"PRIu64" [MB/s]\n", ds / dt);
+#endif
+	}
+#endif
+
 	free(buffer_alloc);
 	free(buffer);
 }
diff --git a/state.c b/state.c
index 494de621..320df9e1 100644
--- a/state.c
+++ b/state.c
@@ -30,6 +30,7 @@ const char* lev_name(unsigned l)
 	case 0 : return "Parity";
 	case 1 : return "Q-Parity";
 	case 2 : return "R-Parity";
+	case 3 : return "S-Parity";
 	}
 
 	return 0;
@@ -41,6 +42,7 @@ const char* lev_config_name(unsigned l)
 	case 0 : return "parity";
 	case 1 : return "q-parity";
 	case 2 : return "r-parity";
+	case 3 : return "s-parity";
 	}
 
 	return 0;
@@ -328,7 +330,7 @@ void state_config(struct snapraid_state* state, const char* path, const char* co
 				exit(EXIT_FAILURE);
 			}
 			state->block_size *= 1024;
-		} else if (strcmp(tag, "parity") == 0 || strcmp(tag, "q-parity") == 0 || strcmp(tag, "r-parity") == 0) {
+		} else if (strcmp(tag, "parity") == 0 || strcmp(tag, "q-parity") == 0 || strcmp(tag, "r-parity") == 0 || strcmp(tag, "s-parity") == 0) {
 			char device[PATH_MAX];
 			char* slash;
 			struct stat st;
@@ -338,6 +340,7 @@ void state_config(struct snapraid_state* state, const char* path, const char* co
 			case 'p' : l = 0; break;
 			case 'q' : l = 1; break;
 			case 'r' : l = 2; break;
+			case 's' : l = 3; break;
 			default:
 				fprintf(stderr, "Invalid '%s' specification in '%s' at line %u\n", tag, path, line);
 				exit(EXIT_FAILURE);
diff --git a/state.h b/state.h
index 6a0bd5e4..871ec87a 100644
--- a/state.h
+++ b/state.h
@@ -24,15 +24,15 @@
 /**
  * Max level of parity supported.
  */
-#define LEV_MAX 3
+#define LEV_MAX 4
 
 /**
- * Returns the parity name: Parity, Q-Parity, R-Parity.
+ * Returns the parity name: Parity, Q-Parity, R-Parity, S-Parity.
  */
 const char* lev_name(unsigned level);
 
 /**
- * Returns the parity name used in the config file: parity, q-parity, r-parity.
+ * Returns the parity name used in the config file: parity, q-parity, r-parity, s-parity.
  */
 const char* lev_config_name(unsigned level);
 
diff --git a/tables.h b/tables.h
index b710e19b..927a8820 100644
--- a/tables.h
+++ b/tables.h
@@ -8795,6 +8795,42 @@ const unsigned char __attribute__((aligned(256))) gfexp4[256] =
 	0x58, 0x7d, 0xe9, 0x83, 0x36, 0xd8, 0x47, 0x01,
 };
 
+const unsigned char __attribute__((aligned(256))) gfexp8[256] =
+{
+	0x01, 0x08, 0x40, 0x3a, 0xcd, 0x26, 0x2d, 0x75,
+	0x8f, 0x0c, 0x60, 0x27, 0x25, 0x35, 0xb5, 0xc1,
+	0x46, 0x0a, 0x50, 0xba, 0xb9, 0xa1, 0x61, 0x2f,
+	0x65, 0x0f, 0x78, 0xe7, 0x6b, 0x7f, 0xdf, 0xb6,
+	0xd9, 0x86, 0x44, 0x1a, 0xd0, 0xce, 0x3e, 0xed,
+	0x3b, 0xc5, 0x66, 0x17, 0xb8, 0xa9, 0x21, 0x15,
+	0xa8, 0x29, 0x55, 0x92, 0xe4, 0x73, 0xbf, 0x91,
+	0xfc, 0xb3, 0xf1, 0xdb, 0x96, 0xc4, 0x6e, 0x57,
+	0x82, 0x64, 0x07, 0x38, 0xdd, 0xa6, 0x59, 0xf2,
+	0xc3, 0x56, 0x8a, 0x24, 0x3d, 0xf5, 0xfb, 0x8b,
+	0x2c, 0x7d, 0xcf, 0x36, 0xad, 0x01, 0x08, 0x40,
+	0x3a, 0xcd, 0x26, 0x2d, 0x75, 0x8f, 0x0c, 0x60,
+	0x27, 0x25, 0x35, 0xb5, 0xc1, 0x46, 0x0a, 0x50,
+	0xba, 0xb9, 0xa1, 0x61, 0x2f, 0x65, 0x0f, 0x78,
+	0xe7, 0x6b, 0x7f, 0xdf, 0xb6, 0xd9, 0x86, 0x44,
+	0x1a, 0xd0, 0xce, 0x3e, 0xed, 0x3b, 0xc5, 0x66,
+	0x17, 0xb8, 0xa9, 0x21, 0x15, 0xa8, 0x29, 0x55,
+	0x92, 0xe4, 0x73, 0xbf, 0x91, 0xfc, 0xb3, 0xf1,
+	0xdb, 0x96, 0xc4, 0x6e, 0x57, 0x82, 0x64, 0x07,
+	0x38, 0xdd, 0xa6, 0x59, 0xf2, 0xc3, 0x56, 0x8a,
+	0x24, 0x3d, 0xf5, 0xfb, 0x8b, 0x2c, 0x7d, 0xcf,
+	0x36, 0xad, 0x01, 0x08, 0x40, 0x3a, 0xcd, 0x26,
+	0x2d, 0x75, 0x8f, 0x0c, 0x60, 0x27, 0x25, 0x35,
+	0xb5, 0xc1, 0x46, 0x0a, 0x50, 0xba, 0xb9, 0xa1,
+	0x61, 0x2f, 0x65, 0x0f, 0x78, 0xe7, 0x6b, 0x7f,
+	0xdf, 0xb6, 0xd9, 0x86, 0x44, 0x1a, 0xd0, 0xce,
+	0x3e, 0xed, 0x3b, 0xc5, 0x66, 0x17, 0xb8, 0xa9,
+	0x21, 0x15, 0xa8, 0x29, 0x55, 0x92, 0xe4, 0x73,
+	0xbf, 0x91, 0xfc, 0xb3, 0xf1, 0xdb, 0x96, 0xc4,
+	0x6e, 0x57, 0x82, 0x64, 0x07, 0x38, 0xdd, 0xa6,
+	0x59, 0xf2, 0xc3, 0x56, 0x8a, 0x24, 0x3d, 0xf5,
+	0xfb, 0x8b, 0x2c, 0x7d, 0xcf, 0x36, 0xad, 0x01,
+};
+
 const unsigned char __attribute__((aligned(256))) gfinv[256] =
 {
 	/* note that the first element is not significative */

From 9af05bce98b985e0f2077a839203dcb49faba207 Mon Sep 17 00:00:00 2001
From: Andrea Mazzoleni <amadvance@gmail.com>
Date: Wed, 18 Sep 2013 21:56:06 +0200
Subject: [PATCH 2/4] Uses a more generic strategy for recoving.

---
 check.c | 485 ++++++++++++++++++++++++++------------------------------
 1 file changed, 224 insertions(+), 261 deletions(-)

diff --git a/check.c b/check.c
index c2ee95d3..545024cb 100644
--- a/check.c
+++ b/check.c
@@ -89,6 +89,106 @@ static int blockcmp(struct snapraid_state* state, int rehash, struct snapraid_bl
 	return 0;
 }
 
+/**
+ * Combinations of 1 of 4 elements.
+ */
+static const struct combo1 {
+	int a;
+} COMBO1[] = {
+	{ 0 },
+	{ 1 },
+	{ 2 },
+	{ 3 },
+};
+
+/**
+ * Combinations of 2 of 4 elements.
+ */
+static const struct combo2 {
+	int a;
+	int b;
+} COMBO2[] = {
+	{ 0, 1 },
+	{ 0, 2 },
+	{ 0, 3 },
+	{ 1, 2 },
+	{ 1, 3 },
+	{ 2, 3 }
+};
+
+/**
+ * Combinations of 3 of 4 elements.
+ */
+static const struct combo3 {
+	int a;
+	int b;
+	int c;
+} COMBO3[] = {
+	{ 0, 1, 2 },
+	{ 0, 1, 3 },
+	{ 0, 2, 3 },
+	{ 1, 2, 3 }
+};
+
+/**
+ * Combinations of 4 of 4 elements.
+ */
+static const struct combo4 {
+	int a;
+	int b;
+	int c;
+	int d;
+} COMBO4[] = {
+	{ 0, 1, 2, 4 },
+};
+
+/**
+ * Checks if the hash of at least one failed block is now matching.
+ */
+static int is_hash_matching(struct snapraid_state* state, int rehash, unsigned diskmax, struct failed_struct* failed, unsigned* failed_map, unsigned failed_count, unsigned char** buffer, unsigned char* buffer_zero)
+{
+	unsigned j;
+	int hash_checked;
+
+	hash_checked = 0; /* keep track if we check at least one block */
+
+	/* check if the recovered blocks are OK */
+	for(j=0;j<failed_count;++j) {
+		if (block_has_updated_hash(failed[failed_map[j]].block)) { /* if the block can be checked */
+			hash_checked = 1;
+			if (blockcmp(state, rehash, failed[failed_map[j]].block, buffer[failed[failed_map[j]].index], buffer_zero) != 0)
+				break;
+		}
+	}
+
+	/* if we checked something, and no block failed the check */
+	if (hash_checked && j==failed_count) {
+		/* recompute all the redundancy information */
+		raid_gen(state->level, buffer, diskmax, state->block_size);
+		return 1;
+	}
+
+	return 0;
+}
+
+/**
+ * Checks if specified parity is now matching with a recomputed one.
+ */
+static int is_parity_matching(struct snapraid_state* state, unsigned diskmax, unsigned i, unsigned char** buffer, unsigned char** buffer_recov)
+{
+	/* recompute parity, note that we don't need parity over i */
+	raid_gen(i + 1, buffer, diskmax, state->block_size);
+
+	/* if the recovered parity block matches */
+	if (memcmp(buffer[diskmax+i], buffer_recov[i], state->block_size) == 0) {
+		/* recompute all the redundancy information */
+		raid_gen(state->level, buffer, diskmax, state->block_size);
+		return 1;
+	}
+
+	return 0;
+}
+
 /**
  * Repair errors.
  * Returns <0 if failure for missing strategy, >0 if data is wrong and we cannot rebuild correctly, 0 on success.
@@ -96,13 +196,9 @@ static int blockcmp(struct snapraid_state* state, int rehash, struct snapraid_bl
  */
 static int repair_step(struct snapraid_state* state, int rehash, unsigned pos, unsigned diskmax, struct failed_struct* failed, unsigned* failed_map, unsigned failed_count, unsigned char** buffer, unsigned char** buffer_recov, unsigned char* buffer_zero)
 {
-	unsigned j;
-	int error = 0;
+	unsigned i, j;
+	int error;
 	int has_hash;
-	int hash_checked;
-	unsigned char* buffer_parity = buffer_recov[0];
-	unsigned char* buffer_qarity = buffer_recov[1];
-	unsigned char* buffer_rarity = buffer_recov[2];
 
 	/* no fix required */
 	if (failed_count == 0) {
@@ -111,6 +207,8 @@ static int repair_step(struct snapraid_state* state, int rehash, unsigned pos, u
 		return 0;
 	}
 
+	error = 0;
+
 	/* check if there is at least a failed block that can be checked for correctness using the hash */
 	/* if there isn't, we have to sacrifice a parity block to check that the result is correct */
 	has_hash = 0;
@@ -119,288 +217,160 @@ static int repair_step(struct snapraid_state* state, int rehash, unsigned pos, u
 			has_hash = 1;
 	}
 
-	/* RAID5 recovering from p */
-	if (failed_count == 1 && buffer_parity != 0 && has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax], buffer_parity, state->block_size);
+	/* recover one block with hash */
+	if (failed_count == 1 && has_hash) {
+		for(i=0;i<sizeof(COMBO1)/sizeof(COMBO1[0]);++i) {
+			/* if parity is missing, do nothing */
+			if (buffer_recov[COMBO1[i].a] == 0)
+				continue;
 
-		/* recover data */
-		raid5_recov_data(buffer, diskmax, state->block_size, failed[failed_map[0]].index);
+			/* copy the redundancy to use */
+			memcpy(buffer[diskmax+COMBO1[i].a], buffer_recov[COMBO1[i].a], state->block_size);
 
-		hash_checked = 0; /* keep track if we check at least one block */
+			/* recover data */
+			raid_recov_1data(buffer, diskmax, state->block_size, failed[failed_map[0]].index, COMBO1[i].a, buffer_zero);
 
-		/* check if the recovered blocks are OK */
-		for(j=0;j<failed_count;++j) {
-			if (block_has_updated_hash(failed[failed_map[j]].block)) { /* if the block can be checked */
-				hash_checked = 1;
-				if (blockcmp(state, rehash, failed[failed_map[j]].block, buffer[failed[failed_map[j]].index], buffer_zero) != 0)
-					break;
-			}
+			if (is_hash_matching(state, rehash, diskmax, failed, failed_map, failed_count, buffer, buffer_zero))
+				return 0;
+
+			fprintf(stdlog, "error:%u:%s: Data error\n", pos, lev_config_name(COMBO1[i].a));
+			++error;
 		}
-
-		/* if we checked something, and no block failed the check */
-		if (hash_checked && j==failed_count) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
-		}
-
-		fprintf(stdlog, "error:%u:parity: Data error\n", pos);
-		++error;
 	}
 
-	/* RAID6 recovering from q, ignoring p */
-	if (failed_count == 1 && buffer_qarity != 0 && has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax+1], buffer_qarity, state->block_size);
+	/* recover one block without hash  */
+	if (failed_count == 1 && !has_hash) {
+		for(i=0;i<sizeof(COMBO2)/sizeof(COMBO2[0]);++i) {
+			/* if parity is missing, do nothing */
+			if (buffer_recov[COMBO2[i].a] == 0 || buffer_recov[COMBO2[i].b] == 0)
+				continue;
 
-		/* recover data */
-		raid6_recov_datap(buffer, diskmax, state->block_size, failed[failed_map[0]].index, buffer_zero);
+			/* copy the redundancy to use */
+			memcpy(buffer[diskmax+COMBO2[i].a], buffer_recov[COMBO2[i].a], state->block_size);
 
-		hash_checked = 0; /* keep track if we check at least one block */
+			/* recover data */
+			raid_recov_1data(buffer, diskmax, state->block_size, failed[failed_map[0]].index, COMBO2[i].a, buffer_zero);
 
-		/* check if the recovered blocks are OK */
-		for(j=0;j<failed_count;++j) {
-			if (block_has_updated_hash(failed[failed_map[j]].block)) { /* if the block can be checked */
-				hash_checked = 1;
-				if (blockcmp(state, rehash, failed[failed_map[j]].block, buffer[failed[failed_map[j]].index], buffer_zero) != 0)
-					break;
-			}
+			if (is_parity_matching(state, diskmax, COMBO2[i].b, buffer, buffer_recov))
+				return 0;
+
+			fprintf(stdlog, "error:%u:%s/%s: Data error\n", pos, lev_config_name(COMBO2[i].a), lev_config_name(COMBO2[i].b));
+			++error;
 		}
-
-		/* if we checked something, and no block failed the check */
-		if (hash_checked && j==failed_count) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
-		}
-
-		fprintf(stdlog, "error:%u:q-parity: Data error\n", pos);
-		++error;
 	}
 
-	/* RAID6 recovering from q, using p to check the result */
-	if (failed_count == 1 && buffer_parity != 0 && buffer_qarity != 0 && !has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax+1], buffer_qarity, state->block_size);
+	/* recover two blocks with hash */
+	if (failed_count == 2 && has_hash) {
+		for(i=0;i<sizeof(COMBO2)/sizeof(COMBO2[0]);++i) {
+			/* if parity is missing, do nothing */
+			if (buffer_recov[COMBO2[i].a] == 0 || buffer_recov[COMBO2[i].b] == 0)
+				continue;
 
-		/* recover data */
-		raid6_recov_datap(buffer, diskmax, state->block_size, failed[failed_map[0]].index, buffer_zero);
+			/* copy the redundancy to use */
+			memcpy(buffer[diskmax+COMBO2[i].a], buffer_recov[COMBO2[i].a], state->block_size);
+			memcpy(buffer[diskmax+COMBO2[i].b], buffer_recov[COMBO2[i].b], state->block_size);
 
-		/* recompute parity with RAID5 to get p from it */
-		raid_gen(1, buffer, diskmax, state->block_size);
+			/* recover data */
+			raid_recov_2data(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, COMBO2[i].a, COMBO2[i].b, buffer_zero);
 
-		/* if the recovered parity block matches */
-		if (memcmp(buffer[diskmax], buffer_parity, state->block_size) == 0) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
+			if (is_hash_matching(state, rehash, diskmax, failed, failed_map, failed_count, buffer, buffer_zero))
+				return 0;
+
+			fprintf(stdlog, "error:%u:%s/%s: Data error\n", pos, lev_config_name(COMBO2[i].a), lev_config_name(COMBO2[i].b));
+			++error;
 		}
-
-		fprintf(stdlog, "error:%u:parity/q-parity: Data error\n", pos);
-		++error;
 	}
 
-	/* RAID6 recovering from p+q */
-	if (failed_count == 2 && buffer_parity != 0 && buffer_qarity != 0 && has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax], buffer_parity, state->block_size);
-		memcpy(buffer[diskmax+1], buffer_qarity, state->block_size);
+	/* recover two blocks without hash */
+	if (failed_count == 2 && !has_hash) {
+		for(i=0;i<sizeof(COMBO3)/sizeof(COMBO3[0]);++i) {
+			/* if parity is missing, do nothing */
+			if (buffer_recov[COMBO3[i].a] == 0 || buffer_recov[COMBO3[i].b] == 0 || buffer_recov[COMBO3[i].c] == 0)
+				continue;
 
-		/* recover 2 data */
-		raid6_recov_2data(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, buffer_zero);
+			/* copy the redundancy to use */
+			memcpy(buffer[diskmax+COMBO3[i].a], buffer_recov[COMBO3[i].a], state->block_size);
+			memcpy(buffer[diskmax+COMBO3[i].b], buffer_recov[COMBO3[i].b], state->block_size);
 
-		hash_checked = 0; /* keep track if we check at least one block */
+			/* recover data */
+			raid_recov_2data(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, COMBO3[i].a, COMBO3[i].b, buffer_zero);
 
-		/* check if the recovered blocks are OK */
-		for(j=0;j<failed_count;++j) {
-			if (block_has_updated_hash(failed[failed_map[j]].block)) { /* if the block can be checked */
-				hash_checked = 1;
-				if (blockcmp(state, rehash, failed[failed_map[j]].block, buffer[failed[failed_map[j]].index], buffer_zero) != 0)
-					break;
-			}
+			if (is_parity_matching(state, diskmax, COMBO3[i].c, buffer, buffer_recov))
+				return 0;
+
+			fprintf(stdlog, "error:%u:%s/%s/%s: Data error\n", pos, lev_config_name(COMBO3[i].a), lev_config_name(COMBO3[i].b), lev_config_name(COMBO3[i].c));
+			++error;
 		}
-
-		/* if we checked something, and no block failed the check */
-		if (hash_checked && j==failed_count) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
-		}
-
-		fprintf(stdlog, "error:%u:parity/q-parity: Data error\n", pos);
-		++error;
 	}
 
-	/* RAIDTP recovering from r, ignoring p,q */
-	if (failed_count == 1 && buffer_rarity != 0 && has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax+2], buffer_rarity, state->block_size);
+	/* recover three blocks with hash */
+	if (failed_count == 3 && has_hash) {
+		for(i=0;i<sizeof(COMBO3)/sizeof(COMBO3[0]);++i) {
+			/* if parity is missing, do nothing */
+			if (buffer_recov[COMBO3[i].a] == 0 || buffer_recov[COMBO3[i].b] == 0 || buffer_recov[COMBO3[i].c] == 0)
+				continue;
 
-		/* recover data */
-		raidTP_recov_datapq(buffer, diskmax, state->block_size, failed[failed_map[0]].index, buffer_zero);
+			/* copy the redundancy to use */
+			memcpy(buffer[diskmax+COMBO3[i].a], buffer_recov[COMBO3[i].a], state->block_size);
+			memcpy(buffer[diskmax+COMBO3[i].b], buffer_recov[COMBO3[i].b], state->block_size);
+			memcpy(buffer[diskmax+COMBO3[i].c], buffer_recov[COMBO3[i].c], state->block_size);
 
-		hash_checked = 0; /* keep track if we check at least one block */
+			/* recover data */
+			raid_recov_3data(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, failed[failed_map[2]].index, COMBO3[i].a, COMBO3[i].b, COMBO3[i].c, buffer_zero);
 
-		/* check if the recovered blocks are OK */
-		for(j=0;j<failed_count;++j) {
-			if (block_has_updated_hash(failed[failed_map[j]].block)) { /* if the block can be checked */
-				hash_checked = 1;
-				if (blockcmp(state, rehash, failed[failed_map[j]].block, buffer[failed[failed_map[j]].index], buffer_zero) != 0)
-					break;
-			}
+			if (is_hash_matching(state, rehash, diskmax, failed, failed_map, failed_count, buffer, buffer_zero))
+				return 0;
+
+			fprintf(stdlog, "error:%u:%s/%s/%s: Data error\n", pos, lev_config_name(COMBO3[i].a), lev_config_name(COMBO3[i].b), lev_config_name(COMBO3[i].c));
+			++error;
 		}
-
-		/* if we checked something, and no block failed the check */
-		if (hash_checked && j==failed_count) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
-		}
-
-		fprintf(stdlog, "error:%u:rarity: Data error\n", pos);
-		++error;
 	}
 
-	/* RAIDTP recovering from r, using p or q to check the result */
-	if (failed_count == 1 && (buffer_parity != 0 || buffer_qarity != 0) && buffer_rarity != 0 && !has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax+2], buffer_rarity, state->block_size);
+	/* recover three blocks without hash */
+	if (failed_count == 3 && !has_hash) {
+		for(i=0;i<sizeof(COMBO4)/sizeof(COMBO4[0]);++i) {
+			/* if parity is missing, do nothing */
+			if (buffer_recov[COMBO4[i].a] == 0 || buffer_recov[COMBO4[i].b] == 0 || buffer_recov[COMBO4[i].c] == 0 || buffer_recov[COMBO4[i].d] == 0)
+				continue;
 
-		/* recover data */
-		raidTP_recov_datapq(buffer, diskmax, state->block_size, failed[failed_map[0]].index, buffer_zero);
+			/* copy the redundancy to use */
+			memcpy(buffer[diskmax+COMBO4[i].a], buffer_recov[COMBO4[i].a], state->block_size);
+			memcpy(buffer[diskmax+COMBO4[i].b], buffer_recov[COMBO4[i].b], state->block_size);
+			memcpy(buffer[diskmax+COMBO4[i].c], buffer_recov[COMBO4[i].c], state->block_size);
 
-		/* recompute parity with RAID6 to get p and q from it */
-		raid_gen(2, buffer, diskmax, state->block_size);
+			/* recover data */
+			raid_recov_3data(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, failed[failed_map[2]].index, COMBO4[i].a, COMBO4[i].b, COMBO4[i].c, buffer_zero);
 
-		/* if the recovered parity block matches */
-		if ((buffer_parity != 0 && memcmp(buffer[diskmax], buffer_parity, state->block_size) == 0)
-			|| (buffer_qarity != 0 && memcmp(buffer[diskmax+1], buffer_qarity, state->block_size) == 0)
-		) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
+			if (is_parity_matching(state, diskmax, COMBO4[i].d, buffer, buffer_recov))
+				return 0;
+
+			fprintf(stdlog, "error:%u:%s/%s/%s/%s: Data error\n", pos, lev_config_name(COMBO4[i].a), lev_config_name(COMBO4[i].b), lev_config_name(COMBO4[i].c), lev_config_name(COMBO4[i].d));
+			++error;
 		}
-
-		fprintf(stdlog, "error:%u:parity/qarity/rarity: Data error\n", pos);
-		++error;
 	}
 
-	/* RAIDTP recovering from p+r */
-	if (failed_count == 2 && buffer_parity != 0 && buffer_rarity != 0 && has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax], buffer_parity, state->block_size);
-		memcpy(buffer[diskmax+2], buffer_rarity, state->block_size);
+	/* recover four blocks with hash */
+	if (failed_count == 4 && has_hash) {
+		for(i=0;i<sizeof(COMBO4)/sizeof(COMBO4[0]);++i) {
+			/* if parity is missing, do nothing */
+			if (buffer_recov[COMBO4[i].a] == 0 || buffer_recov[COMBO4[i].b] == 0 || buffer_recov[COMBO4[i].c] == 0 || buffer_recov[COMBO4[i].d] == 0)
+				continue;
 
-		/* recover 2 data */
-		raidTP_recov_2dataq(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, buffer_zero);
+			/* copy the redundancy to use */
+			memcpy(buffer[diskmax+COMBO4[i].a], buffer_recov[COMBO4[i].a], state->block_size);
+			memcpy(buffer[diskmax+COMBO4[i].b], buffer_recov[COMBO4[i].b], state->block_size);
+			memcpy(buffer[diskmax+COMBO4[i].c], buffer_recov[COMBO4[i].c], state->block_size);
+			memcpy(buffer[diskmax+COMBO4[i].d], buffer_recov[COMBO4[i].d], state->block_size);
 
-		hash_checked = 0; /* keep track if we check at least one block */
+			/* recover data */
+			raid_recov_4data(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, failed[failed_map[2]].index, failed[failed_map[3]].index, COMBO4[i].a, COMBO4[i].b, COMBO4[i].c, COMBO4[i].d, buffer_zero);
 
-		/* check if the recovered blocks are OK */
-		for(j=0;j<failed_count;++j) {
-			if (block_has_updated_hash(failed[failed_map[j]].block)) { /* if the block can be checked */
-				hash_checked = 1;
-				if (blockcmp(state, rehash, failed[failed_map[j]].block, buffer[failed[failed_map[j]].index], buffer_zero) != 0)
-					break;
-			}
+			if (is_hash_matching(state, rehash, diskmax, failed, failed_map, failed_count, buffer, buffer_zero))
+				return 0;
+
+			fprintf(stdlog, "error:%u:%s/%s/%s/%s: Data error\n", pos, lev_config_name(COMBO4[i].a), lev_config_name(COMBO4[i].b), lev_config_name(COMBO4[i].c), lev_config_name(COMBO4[i].d));
+			++error;
 		}
-
-		/* if we checked something, and no block failed the check */
-		if (hash_checked && j==failed_count) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
-		}
-
-		fprintf(stdlog, "error:%u:parity/rarity: Data error\n", pos);
-		++error;
-	}
-
-	/* RAIDTP recovering from q+r */
-	if (failed_count == 2 && buffer_qarity != 0 && buffer_rarity != 0 && has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax+1], buffer_qarity, state->block_size);
-		memcpy(buffer[diskmax+2], buffer_rarity, state->block_size);
-
-		/* recover 2 data */
-		raidTP_recov_2datap(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, buffer_zero);
-
-		hash_checked = 0; /* keep track if we check at least one block */
-
-		/* check if the recovered blocks are OK */
-		for(j=0;j<failed_count;++j) {
-			if (block_has_updated_hash(failed[failed_map[j]].block)) { /* if the block can be checked */
-				hash_checked = 1;
-				if (blockcmp(state, rehash, failed[failed_map[j]].block, buffer[failed[failed_map[j]].index], buffer_zero) != 0)
-					break;
-			}
-		}
-
-		/* if we checked something, and no block failed the check */
-		if (hash_checked && j==failed_count) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
-		}
-
-		fprintf(stdlog, "error:%u:qarity/rarity: Data error\n", pos);
-		++error;
-	}
-
-	/* RAIDTP recovering from q+r, using p to check the result */
-	if (failed_count == 2 && buffer_parity != 0 && buffer_qarity != 0 && buffer_rarity != 0 && !has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax+1], buffer_qarity, state->block_size);
-		memcpy(buffer[diskmax+2], buffer_rarity, state->block_size);
-
-		/* recover data */
-		raidTP_recov_2datap(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, buffer_zero);
-
-		/* recompute parity with RAID5 to get p from it */
-		raid_gen(1, buffer, diskmax, state->block_size);
-
-		/* if the recovered parity block matches */
-		if (buffer_parity != 0 && memcmp(buffer[diskmax], buffer_parity, state->block_size) == 0) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
-		}
-
-		fprintf(stdlog, "error:%u:parity/qarity/rarity: Data error\n", pos);
-		++error;
-	}
-
-	/* RAIDTP recovering from p+q+r */
-	if (failed_count == 3 && buffer_parity != 0 && buffer_qarity != 0 && buffer_rarity != 0 && has_hash) {
-		/* copy the redundancy to use */
-		memcpy(buffer[diskmax], buffer_parity, state->block_size);
-		memcpy(buffer[diskmax+1], buffer_qarity, state->block_size);
-		memcpy(buffer[diskmax+2], buffer_rarity, state->block_size);
-
-		/* recover 3 data */
-		raidTP_recov_3data(buffer, diskmax, state->block_size, failed[failed_map[0]].index, failed[failed_map[1]].index, failed[failed_map[2]].index, buffer_zero);
-
-		hash_checked = 0; /* keep track if we check at least one block */
-
-		/* check if the recovered blocks are OK */
-		for(j=0;j<failed_count;++j) {
-			if (block_has_updated_hash(failed[failed_map[j]].block)) { /* if the block can be checked */
-				hash_checked = 1;
-				if (blockcmp(state, rehash, failed[failed_map[j]].block, buffer[failed[failed_map[j]].index], buffer_zero) != 0)
-					break;
-			}
-		}
-
-		/* if we checked something, and no block failed the check */
-		if (hash_checked && j==failed_count) {
-			/* recompute all the redundancy information */
-			raid_gen(state->level, buffer, diskmax, state->block_size);
-			return 0;
-		}
-
-		fprintf(stdlog, "error:%u:parity/qarity/rarity: Data error\n", pos);
-		++error;
 	}
 
 	/* return the number of failed attempts, or -1 if no strategy */
@@ -914,19 +884,12 @@ static int state_check_process(struct snapraid_state* state, int check, int fix,
 			unsigned char* buffer_zero;
 
 			/* buffers for parity read and not computed */
-			if (state->level == 1) {
-				buffer_recov[0] = buffer[diskmax + 1];
-				buffer_recov[1] = 0;
-				buffer_recov[2] = 0;
-			} else if (state->level == 2) {
-				buffer_recov[0] = buffer[diskmax + 2];
-				buffer_recov[1] = buffer[diskmax + 3];
-				buffer_recov[2] = 0;
-			} else {
-				buffer_recov[0] = buffer[diskmax + 3];
-				buffer_recov[1] = buffer[diskmax + 4];
-				buffer_recov[2] = buffer[diskmax + 5];
-			}
+			for(l=0;l<state->level;++l)
+				buffer_recov[l] = buffer[diskmax + state->level + l];
+			for(;l<LEV_MAX;++l)
+				buffer_recov[l] = 0;
+
+			/* the zero buffer is the last one */
 			buffer_zero = buffer[buffermax-1];
 
 			/* read the parity */

From f4c1e288b2d52ee38df5ed5f10ae5d34c462a655 Mon Sep 17 00:00:00 2001
From: Andrea Mazzoleni <amadvance@gmail.com>
Date: Wed, 18 Sep 2013 23:06:07 +0200
Subject: [PATCH 3/4] Removes the specific recovering functions. We now always
 use the generic ones.

---
 raid.c     | 261 ++++++++++++++-----------------------
 raid.h     |  44 +------
 selftest.c | 368 +++--------------------------------------------------
 3 files changed, 122 insertions(+), 551 deletions(-)

diff --git a/raid.c b/raid.c
index bd470a85..e584d61e 100644
--- a/raid.c
+++ b/raid.c
@@ -1232,8 +1232,8 @@ void raid_gen(unsigned level, unsigned char** buffer, unsigned diskmax, unsigned
  * This method is also the same used by ZFS to implement its RAIDTP
  * support.
  *
- * Note that the same extension is not possible for Quad Parity because
- * in the approach used in the paper we don't have the guarantee to always have
+ * To support RAIDQP (Quad Parity) we go futher using also the parity generator "8".
+ * Note that in this case we don't have the guarantee to always have
  * a system of independent linear equations, and for Quad Parity in some cases
  * the system is not solvable.
  *
@@ -1258,18 +1258,18 @@ void raid_gen(unsigned level, unsigned char** buffer, unsigned diskmax, unsigned
  * Q = sum(2^i * Di) 0<=i<n
  * R = sum(4^i * Di) 0<=i<n
  *
- * To recover from a failure of three disks at indes x,y,z,
+ * To recover from a failure of three disks at indexes x,y,z,
  * with 0<=x<y<z<n, we compute the syndromes of the available disks:
  *
- * Pxyz = sum(Di) 0<=i<n,i!=x,i!=y,i!=z
- * Qxyz = sum(2^i * Di) 0<=i<n,i!=x,i!=y,i!=z
- * Rxyz = sum(4^i * Di) 0<=i<n,i!=x,i!=y,i!=z
+ * Pa = sum(Di) 0<=i<n,i!=x,i!=y,i!=z
+ * Qa = sum(2^i * Di) 0<=i<n,i!=x,i!=y,i!=z
+ * Ra = sum(4^i * Di) 0<=i<n,i!=x,i!=y,i!=z
  *
  * and if we define:
  *
- * Pd = Pxyz + P
- * Qd = Qxyz + Q
- * Rd = Rxyz + R
+ * Pd = Pa + P
+ * Qd = Qa + Q
+ * Rd = Ra + R
  *
  * we can sum these two set equations, obtaining:
  *
@@ -1277,27 +1277,31 @@ void raid_gen(unsigned level, unsigned char** buffer, unsigned diskmax, unsigned
  * Qd = 2^x * Dx + 2^y * Dy + 2^z * Dz
  * Rd = 4^x * Dx + 4^y * Dy + 4^z * Dz
  *
- * A linear system of three equations solvable by substitution on Dx, Dy and Dz,
- * because all the others variables are known and the equations are independent
- * for any 0<=x<y<z<255.
- * We can prove that the equations are independent by brute-force, trying all the
- * possible combinations of x,y,z.
+ * A linear system always solvable because the coefficients matrix is always
+ * not singular, including all its submatrix.
+ * We can prove that by brute-force, trying all the possible combinations
+ * of x,y,z with 0<=x<y<z<255.
  *
  * The other recovering cases are a simplification of the general one,
  * with some equations or addends removed.
  *
- * We have in total 7 different recovering strategies:
+ * For Quad Parity we follow the same method starting with:
  *
- * RAID5 - Dx recovered with P
- * RAID6 - Dx recovered with Q
- * RAID6 - Dx and Dy recovered with P,Q
- * RAIDTP - Dx recovered with R
- * RAIDTP - Dx,Dy recovered with P,R
- * RAIDTP - Dx,Dy recovered with Q,R
- * RAIDTP - Dx,Dy,Dz recovered with P,Q,R
+ * P = sum(Di) 0<=i<n
+ * Q = sum(2^i * Di) 0<=i<n
+ * R = sum(4^i * Di) 0<=i<n
+ * S = sum(8^i * Di) 0<=i<n
  *
- * each one implemented in a different function using linear system
- * of equations with a different solution.
+ * obtaining:
+ *
+ * Pd =       Dx +       Dy +       Dz +       Dv
+ * Qd = 2^x * Dx + 2^y * Dy + 2^z * Dz + 2^v * Dv
+ * Rd = 4^x * Dx + 4^y * Dy + 4^z * Dz + 4^v * Dv
+ * Sd = 8^x * Dx + 8^y * Dy + 8^z * Dz + 8^v * Dv
+ *
+ * In this case the coefficients matrix and all its submatrix are not singular
+ * only with up to 21 data disk (0<=x<y<z<v<21 with GF primitive poly 285).
+ * Note that the cases failing with 22 data disk are some matrix 3x3.
  *
  * References:
  * [1] Anvin, "The mathematics of RAID-6", 2004
@@ -1390,7 +1394,9 @@ static inline const unsigned char* table(unsigned char a)
 	return gfmul[a];
 }
 
-/*
+/**
+ * Recover failure of one data block for RAID5.
+ *
  * Starting from the equation:
  *
  * Pd = Dx
@@ -1399,7 +1405,7 @@ static inline const unsigned char* table(unsigned char a)
  *
  * Dx = Pd (this one is easy :D)
  */
-void raid5_recov_data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x)
+static void raid5_recov_data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x)
 {
 	unsigned char* p;
 	unsigned char* dp;
@@ -1418,7 +1424,9 @@ void raid5_recov_data(unsigned char** dptrs, unsigned diskmax, unsigned size, in
 	dptrs[diskmax] = p;
 }
 
-/*
+/**
+ * Recover failure of two data blocks for RAID6.
+ *
  * Starting from the equations:
  *
  * Pd = Dx + Dy
@@ -1439,7 +1447,7 @@ void raid5_recov_data(unsigned char** dptrs, unsigned diskmax, unsigned size, in
  *
  * That are always satisfied for any 0<=x<y<255.
  */
-void raid6_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, unsigned char* zero)
+static void raid6_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, unsigned char* zero)
 {
 	unsigned char* p;
 	unsigned char* dp;
@@ -1498,120 +1506,9 @@ void raid6_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, i
 	}
 }
 
-/*
- * Starting from the equations:
+/**
+ * Recover failure of two data blocks without Q for RAIDTP.
  *
- * Qd = 2^x * Dx
- *
- * and solving we get:
- *
- * Dx = 2^(-x) * Qd
- *
- * with conditions:
- *
- * 2^x != 0
- *
- * That are always satisfied for any 0<=x<y<255.
- */
-void raid6_recov_datap(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, unsigned char* zero)
-{
-	unsigned char* q;
-	unsigned char* dq;
-	const unsigned char* qxmul; /* Q multiplier to compute Dx */
-
-	q = dptrs[diskmax+1];
-
-	/* compute syndrome with zero for the missing data page */
-	/* use the dead data page as temporary storage for delta q */
-	dq = dptrs[x];
-	dptrs[x] = zero;
-	dptrs[diskmax+1] = dq;
-
-	raid6_gen(dptrs, diskmax, size);
-
-	/* restore pointers */
-	dptrs[x] = dq;
-	dptrs[diskmax+1] = q;
-
-	/* select tables */
-	qxmul = table( inv(pow2(x)) );
-
-	while (size--) {
-		/* delta */
-		unsigned char Qd = *q ^ *dq;
-
-		/* addends to reconstruct Dx */
-		unsigned char qam = qxmul[Qd];
-
-		/* reconstruct Dx */
-		unsigned char Dx = qam;
-
-		/* set */
-		*dq = Dx;
-
-		++q;
-		++dq;
-	}
-}
-
-/*
- * Starting from the equation:
- *
- * Rd = 4^x * Dx
- *
- * and solving we get:
- *
- * Dx = 4^(-x) * Rd
- *
- * with conditions:
- *
- * 4^x != 0
- *
- * That are always satisfied for any 0<=x<y<255.
- *
- */
-void raidTP_recov_datapq(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, unsigned char* zero)
-{
-	unsigned char* r;
-	unsigned char* dr;
-	const unsigned char* rxmul; /* R multiplier to compute Dx */
-
-	r = dptrs[diskmax+2];
-
-	/* compute syndrome with zero for the missing data page */
-	/* use the dead data page as temporary storage for delta r */
-	dr = dptrs[x];
-	dptrs[x] = zero;
-	dptrs[diskmax+2] = dr;
-
-	raidTP_gen(dptrs, diskmax, size);
-
-	/* restore pointers */
-	dptrs[x] = dr;
-	dptrs[diskmax+2] = r;
-
-	/* select tables */
-	rxmul = table( inv(pow4(x)) );
-
-	while (size--) {
-		/* delta */
-		unsigned char Rd = *r ^ *dr;
-
-		/* addends to reconstruct Dx */
-		unsigned char rxm = rxmul[Rd];
-
-		/* reconstruct Dx */
-		unsigned char Dx = rxm;
-
-		/* set */
-		*dr = Dx;
-
-		++r;
-		++dr;
-	}
-}
-
-/*
  * Starting from the equations:
  *
  * Pd = Dx + Dy
@@ -1632,7 +1529,7 @@ void raidTP_recov_datapq(unsigned char** dptrs, unsigned diskmax, unsigned size,
  *
  * That are always satisfied for any 0<=x<y<255.
  */
-void raidTP_recov_2dataq(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, unsigned char* zero)
+static void raidTP_recov_2dataq(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, unsigned char* zero)
 {
 	unsigned char* p;
 	unsigned char* dp;
@@ -1693,7 +1590,9 @@ void raidTP_recov_2dataq(unsigned char** dptrs, unsigned diskmax, unsigned size,
 	}
 }
 
-/*
+/**
+ * Recover failure of two data blocks without P for RAIDTP.
+ *
  * Starting from the equations:
  *
  * Qd = 2^x * Dx + 2^y * Dy
@@ -1715,7 +1614,7 @@ void raidTP_recov_2dataq(unsigned char** dptrs, unsigned diskmax, unsigned size,
  *
  * That are always satisfied for any 0<=x<y<255.
  */
-void raidTP_recov_2datap(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, unsigned char* zero)
+static void raidTP_recov_2datap(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, unsigned char* zero)
 {
 	unsigned char* q;
 	unsigned char* dq;
@@ -1784,7 +1683,9 @@ void raidTP_recov_2datap(unsigned char** dptrs, unsigned diskmax, unsigned size,
 	}
 }
 
-/*
+/**
+ * Recover failure of three data blocks for RAIDTP.
+ *
  * Starting from the equations:
  *
  * Pd = Dx + Dy + Dz
@@ -1817,7 +1718,7 @@ void raidTP_recov_2datap(unsigned char** dptrs, unsigned diskmax, unsigned size,
  *
  * That are always satisfied for any 0<=x<y<z<255.
  */
-void raidTP_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, int z, unsigned char* zero)
+static void raidTP_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size, int x, int y, int z, unsigned char* zero)
 {
 	unsigned char* p;
 	unsigned char* dp;
@@ -1972,18 +1873,20 @@ static void invert(unsigned char* M, unsigned char* V, unsigned n)
 	}
 }
 
-/*
+/**
+ * Recover failure of one data block x using parity i for any RAID level.
+ *
  * Starting from the equation:
  *
- * Sd = 2^i^x * Dx
+ * Sd = (2^i)^x * Dx
  *
  * and solving we get:
  *
- * Dx = 2^i^(-x) * Sd
+ * Dx = (2^i)^(-x) * Sd
  *
  * with conditions:
  *
- * 2^i^x != 0
+ * (2^i)^x != 0
  *
  * That are always satisfied for any 0<=x<y<255.
  */
@@ -1993,9 +1896,11 @@ void raid_recov_1data(unsigned char** dptrs, unsigned diskmax, unsigned size, in
 	unsigned char* dp;
 	const unsigned char* pxmul; /* P multiplier to compute Dx */
 
-	/* if i is P use the fast one */
-	if (i == 0)
+	/* if it's RAID5 uses the dedicated and faster function */
+	if (i == 0) {
 		raid5_recov_data(dptrs, diskmax, size, x);
+		return;
+	}
 
 	p = dptrs[diskmax+i];
 
@@ -2032,11 +1937,13 @@ void raid_recov_1data(unsigned char** dptrs, unsigned diskmax, unsigned size, in
 	}
 }
 
-/*
+/**
+ * Recover failure of two data blocks x,y using parity i,j for any RAID level.
+ *
  * Starting from the equations:
  *
- * Pd = 2^i^x * Dx + 2^i^y * Dy
- * Qd = 2^j^x * Dx + 2^j^y * Dy
+ * Pd = (2^i)^x * Dx + (2^i)^y * Dy
+ * Qd = (2^j)^x * Dx + (2^j)^y * Dy
  *
  * we solve inverting the coefficients matrix.
  */
@@ -2054,6 +1961,22 @@ void raid_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, in
 	unsigned char V[2*2];
 	unsigned n = 2;
 
+	/* if it's RAID6 uses the dedicated and faster function */
+	if (i == 0 && j == 1) {
+		raid6_recov_2data(dptrs, diskmax, size, x, y, zero);
+		return;
+	}
+	/* if it's RAIDTP uses the dedicated and faster function */
+	if (i == 0 && j == 2) {
+		raidTP_recov_2dataq(dptrs, diskmax, size, x, y, zero);
+		return;
+	}
+	/* if it's RAIDTP uses the dedicated and faster function */
+	if (i == 1 && j == 2) {
+		raidTP_recov_2datap(dptrs, diskmax, size, x, y, zero);
+		return;
+	}
+
 	p = dptrs[diskmax+i];
 	q = dptrs[diskmax+j];
 
@@ -2127,12 +2050,14 @@ void raid_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, in
 	}
 }
 
-/*
+/**
+ * Recover failure of three data blocks x,y,z using parity i,j,k for any RAID level.
+ *
  * Starting from the equations:
  *
- * Pd = 2^i^x * Dx + 2^i^y * Dy + 2^i^z * Dz
- * Qd = 2^j^x * Dx + 2^j^y * Dy + 2^j^z * Dz
- * Rd = 2^k^x * Dx + 2^k^y * Dy + 2^k^z * Dz
+ * Pd = (2^i)^x * Dx + (2^i)^y * Dy + (2^i)^z * Dz
+ * Qd = (2^j)^x * Dx + (2^j)^y * Dy + (2^j)^z * Dz
+ * Rd = (2^k)^x * Dx + (2^k)^y * Dy + (2^k)^z * Dz
  *
  * we solve inverting the coefficients matrix.
  */
@@ -2157,6 +2082,12 @@ void raid_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size, in
 	unsigned char V[3*3];
 	unsigned n = 3;
 
+	/* if it's RAIDTP uses the dedicated and faster function */
+	if (i == 0 && j == 1 && k == 2) {
+		raidTP_recov_3data(dptrs, diskmax, size, x, y, z, zero);
+		return;
+	}
+
 	p = dptrs[diskmax+i];
 	q = dptrs[diskmax+j];
 	r = dptrs[diskmax+k];
@@ -2260,13 +2191,15 @@ void raid_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size, in
 	}
 }
 
-/*
+/**
+ * Recover failure of four data blocks x,y,z,v using parity i,j,k,l for any RAID level.
+ *
  * Starting from the equations:
  *
- * Pd = 2^i^x * Dx + 2^i^y * Dy + 2^i^z * Dz + 2^i^v * Dv
- * Qd = 2^j^x * Dx + 2^j^y * Dy + 2^j^z * Dz + 2^j^v * Dv
- * Rd = 2^k^x * Dx + 2^k^y * Dy + 2^k^z * Dz + 2^k^v * Dv
- * Sd = 2^l^x * Dx + 2^l^y * Dy + 2^l^z * Dz + 2^l^v * Dv
+ * Pd = (2^i)^x * Dx + (2^i)^y * Dy + (2^i)^z * Dz + (2^i)^v * Dv
+ * Qd = (2^j)^x * Dx + (2^j)^y * Dy + (2^j)^z * Dz + (2^j)^v * Dv
+ * Rd = (2^k)^x * Dx + (2^k)^y * Dy + (2^k)^z * Dz + (2^k)^v * Dv
+ * Sd = (2^l)^x * Dx + (2^l)^y * Dy + (2^l)^z * Dz + (2^l)^v * Dv
  *
  * we solve inverting the coefficients matrix.
  */
diff --git a/raid.h b/raid.h
index 120efe68..14fb1194 100644
--- a/raid.h
+++ b/raid.h
@@ -21,50 +21,16 @@
 /****************************************************************************/
 /* raid */
 
+/**
+ * Max number of data disks for quad parity.
+ */
+#define RAID_QUAD_DATA_LIMIT 21
+
 /**
  * Syndrome computation.
  */
 void raid_gen(unsigned level, unsigned char** buffer, unsigned diskmax, unsigned size);
 
-/*
- * Recover failure of one data block.
- */
-void raid5_recov_data(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila);
-
-/*
- * Recover two failed data blocks.
- */
-void raid6_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, unsigned char* zero);
-
-/*
- * Recover failure of one data block ignoring P.
- * Note that P content is not recovered.
- */
-void raid6_recov_datap(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, unsigned char* zero);
-
-/*
- * Recover failure of one data block ignoring P and Q.
- * Note that P and Q content is not recovered.
- */
-void raidTP_recov_datapq(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, unsigned char* zero);
-
-/*
- * Recover failure of two data blocks ignoring P.
- * Note that P content is not recovered.
- */
-void raidTP_recov_2datap(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, unsigned char* zero);
-
-/*
- * Recover failure of two data blocks ignoring Q.
- * Note that Q content is not recovered.
- */
-void raidTP_recov_2dataq(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, unsigned char* zero);
-
-/*
- * Recover failure of three data blocks.
- */
-void raidTP_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, int failc, unsigned char* zero);
-
 /*
  * Recover failure of one data block x using parity i.
  */
diff --git a/selftest.c b/selftest.c
index 01dd4ab3..3e2204b3 100644
--- a/selftest.c
+++ b/selftest.c
@@ -20,6 +20,8 @@
 #include "snapraid.h"
 #include "util.h"
 #include "raid.h"
+#include "elem.h"
+#include "state.h"
 
 struct hash_test_vector {
 	const char* data;
@@ -43,338 +45,7 @@ static struct hash_test_vector TEST_SPOOKY2[] = {
 { 0, 0, { 0 } }
 };
 
-static void raid5test(unsigned diskmax, unsigned block_size)
-{
-	void* buffer_alloc;
-	unsigned char** buffer;
-	unsigned char* buffer_test;
-	unsigned buffermax;
-	unsigned i, j;
-
-	buffermax = diskmax + 1 + 1;
-
-	buffer = malloc_nofail_vector_align(buffermax, block_size, &buffer_alloc);
-
-	buffer_test = buffer[buffermax-1];
-
-	/* fill with random */
-	for(i=0;i<diskmax;++i) {
-		for(j=0;j<block_size;++j)
-			buffer[i][j] = rand();
-	}
-
-	/* compute the parity */
-	raid_gen(1, buffer, diskmax, block_size);
-
-	for(i=0;i<diskmax;++i) {
-		unsigned char* a;
-
-		/* save */
-		a = buffer[i];
-
-		/* damage */
-		buffer[i] = buffer_test;
-
-		/* recover */
-		raid5_recov_data(buffer, diskmax, block_size, i);
-
-		/* restore */
-		buffer[i] = a;
-
-		/* check */
-		if (memcmp(buffer_test, buffer[i], block_size) != 0) {
-			fprintf(stderr, "Failed RAID5 test\n");
-			exit(EXIT_FAILURE);
-		}
-	}
-
-	free(buffer_alloc);
-	free(buffer);
-}
-
-static void raid6test(unsigned diskmax, unsigned block_size)
-{
-	void* buffer_alloc;
-	unsigned char** buffer;
-	unsigned char* buffer_test0;
-	unsigned char* buffer_test1;
-	unsigned char* buffer_zero;
-	unsigned char* buffer_bad;
-	unsigned buffermax;
-	unsigned i, j;
-
-	buffermax = diskmax + 2 + 2 + 1 + 1;
-
-	buffer = malloc_nofail_vector_align(buffermax, block_size, &buffer_alloc);
-
-	buffer_test0 = buffer[buffermax-4];
-	buffer_test1 = buffer[buffermax-3];
-	buffer_zero = buffer[buffermax-2];
-	buffer_bad = buffer[buffermax-1];
-	memset(buffer_zero, 0, block_size);
-	memset(buffer_bad, 0xA5, block_size);
-
-	/* fill with random */
-	for(i=0;i<diskmax;++i) {
-		for(j=0;j<block_size;++j)
-			buffer[i][j] = rand();
-	}
-
-	/* compute the parity */
-	raid_gen(2, buffer, diskmax, block_size);
-
-	/* check 2data */
-	for(i=0;i<diskmax;++i) {
-		for(j=i+1;j<diskmax;++j) {
-			unsigned char* a;
-			unsigned char* b;
-
-			/* save */
-			a = buffer[i];
-			b = buffer[j];
-
-			/* damage */
-			buffer[i] = buffer_test0;
-			buffer[j] = buffer_test1;
-
-			/* recover */
-			raid6_recov_2data(buffer, diskmax, block_size, i, j, buffer_zero);
-
-			/* restore */
-			buffer[i] = a;
-			buffer[j] = b;
-
-			/* check */
-			if (memcmp(buffer_test0, buffer[i], block_size) != 0
-				|| memcmp(buffer_test1, buffer[j], block_size) != 0
-			) {
-				fprintf(stderr, "Failed RAID6 test\n");
-				exit(EXIT_FAILURE);
-			}
-		}
-	}
-
-	/* check datap */
-	for(i=0;i<diskmax;++i) {
-		unsigned char* a;
-		unsigned char* p;
-
-		/* other is p */
-		j = diskmax;
-
-		/* save */
-		a = buffer[i];
-		p = buffer[j];
-
-		/* damage */
-		buffer[i] = buffer_test0;
-		buffer[j] = buffer_bad;
-
-		/* recover */
-		raid6_recov_datap(buffer, diskmax, block_size, i, buffer_zero);
-
-		/* restore */
-		buffer[i] = a;
-		buffer[j] = p;
-
-		/* check */
-		if (memcmp(buffer_test0, buffer[i], block_size) != 0) {
-			fprintf(stderr, "Failed RAID6 test\n");
-			exit(EXIT_FAILURE);
-		}
-	}
-
-	free(buffer_alloc);
-	free(buffer);
-}
-
-static void raidTPtest(unsigned diskmax, unsigned block_size)
-{
-	void* buffer_alloc;
-	unsigned char** buffer;
-	unsigned char* buffer_test0;
-	unsigned char* buffer_test1;
-	unsigned char* buffer_test2;
-	unsigned char* buffer_zero;
-	unsigned char* buffer_bad;
-	unsigned buffermax;
-	unsigned i, j, k;
-
-	buffermax = diskmax + 3 + 3 + 1 + 1;
-
-	buffer = malloc_nofail_vector_align(buffermax, block_size, &buffer_alloc);
-
-	buffer_test0 = buffer[buffermax-5];
-	buffer_test1 = buffer[buffermax-4];
-	buffer_test2 = buffer[buffermax-3];
-	buffer_zero = buffer[buffermax-2];
-	buffer_bad = buffer[buffermax-1];
-	memset(buffer_zero, 0, block_size);
-	memset(buffer_bad, 0xA5, block_size);
-
-	/* fill with random */
-	for(i=0;i<diskmax;++i) {
-		for(j=0;j<block_size;++j)
-			buffer[i][j] = rand();
-	}
-
-	/* compute the parity */
-	raid_gen(3, buffer, diskmax, block_size);
-
-	/* check 3data */
-	for(i=0;i<diskmax;++i) {
-		for(j=i+1;j<diskmax;++j) {
-			for(k=j+1;k<diskmax;++k) {
-				unsigned char* a;
-				unsigned char* b;
-				unsigned char* c;
-
-				/* save */
-				a = buffer[i];
-				b = buffer[j];
-				c = buffer[k];
-
-				/* damage */
-				buffer[i] = buffer_test0;
-				buffer[j] = buffer_test1;
-				buffer[k] = buffer_test2;
-
-				/* recover */
-				raidTP_recov_3data(buffer, diskmax, block_size, i, j, k, buffer_zero);
-
-				/* restore */
-				buffer[i] = a;
-				buffer[j] = b;
-				buffer[k] = c;
-
-				/* check */
-				if (memcmp(buffer_test0, buffer[i], block_size) != 0
-					|| memcmp(buffer_test1, buffer[j], block_size) != 0
-					|| memcmp(buffer_test2, buffer[k], block_size) != 0
-				) {
-					fprintf(stderr, "Failed RAIDTP test\n");
-					exit(EXIT_FAILURE);
-				}
-			}
-		}
-	}
-
-	/* check 2dataq */
-	for(i=0;i<diskmax;++i) {
-		for(j=i+1;j<diskmax;++j) {
-			unsigned char* a;
-			unsigned char* b;
-			unsigned char* q;
-
-			/* other is q */
-			k = diskmax + 1;
-
-			/* save */
-			a = buffer[i];
-			b = buffer[j];
-			q = buffer[k];
-
-			/* damage */
-			buffer[i] = buffer_test0;
-			buffer[j] = buffer_test1;
-			buffer[k] = buffer_bad;
-
-			/* recover */
-			raidTP_recov_2dataq(buffer, diskmax, block_size, i, j, buffer_zero);
-
-			/* restore */
-			buffer[i] = a;
-			buffer[j] = b;
-			buffer[k] = q;
-
-			/* check */
-			if (memcmp(buffer_test0, buffer[i], block_size) != 0
-				|| memcmp(buffer_test1, buffer[j], block_size) != 0
-			) {
-				fprintf(stderr, "Failed RAIDTP test\n");
-				exit(EXIT_FAILURE);
-			}
-		}
-	}
-
-	/* check 2datap */
-	for(i=0;i<diskmax;++i) {
-		for(j=i+1;j<diskmax;++j) {
-			unsigned char* a;
-			unsigned char* b;
-			unsigned char* p;
-
-			/* other is p */
-			k = diskmax;
-
-			/* save */
-			a = buffer[i];
-			b = buffer[j];
-			p = buffer[k];
-
-			/* damage */
-			buffer[i] = buffer_test0;
-			buffer[j] = buffer_test1;
-			buffer[k] = buffer_bad;
-
-			/* recover */
-			raidTP_recov_2datap(buffer, diskmax, block_size, i, j, buffer_zero);
-
-			/* restore */
-			buffer[i] = a;
-			buffer[j] = b;
-			buffer[k] = p;
-
-			/* check */
-			if (memcmp(buffer_test0, buffer[i], block_size) != 0
-				|| memcmp(buffer_test1, buffer[j], block_size) != 0
-			) {
-				fprintf(stderr, "Failed RAIDTP test\n");
-				exit(EXIT_FAILURE);
-			}
-		}
-	}
-
-	/* check datapq */
-	for(i=0;i<diskmax;++i) {
-		unsigned char* a;
-		unsigned char* p;
-		unsigned char* q;
-
-		/* others are p,q */
-		j = diskmax;
-		k = diskmax + 1;
-
-		/* save */
-		a = buffer[i];
-		p = buffer[j];
-		q = buffer[k];
-
-		/* damage */
-		buffer[i] = buffer_test0;
-		buffer[j] = buffer_bad;
-		buffer[k] = buffer_bad;
-
-		/* recover */
-		raidTP_recov_datapq(buffer, diskmax, block_size, i, buffer_zero);
-
-		/* restore */
-		buffer[i] = a;
-		buffer[j] = p;
-		buffer[k] = q;
-
-		/* check */
-		if (memcmp(buffer_test0, buffer[i], block_size) != 0) {
-			fprintf(stderr, "Failed RAIDTP test\n");
-			exit(EXIT_FAILURE);
-		}
-	}
-
-	free(buffer_alloc);
-	free(buffer);
-}
-
-static void raidQPtest(unsigned diskmax, unsigned block_size)
+static void raidtest(unsigned diskmax, unsigned block_size)
 {
 	void* buffer_alloc;
 	unsigned char** buffer;
@@ -391,7 +62,7 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 	struct combo1 {
 		int ignore[3];
 		int use;
-	} COMBO1[4] = {
+	} COMBO1[] = {
 		{ { 1, 2, 3 }, 0 },
 		{ { 0, 2, 3 }, 1 },
 		{ { 0, 1, 3 }, 2 },
@@ -401,7 +72,7 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 	struct combo2 {
 		int ignore[2];
 		int use[2];
-	} COMBO2[6] = {
+	} COMBO2[] = {
 		{ { 0, 1 }, { 2, 3 } },
 		{ { 0, 2 }, { 1, 3 } },
 		{ { 0, 3 }, { 1, 2 } },
@@ -413,14 +84,14 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 	struct combo3 {
 		int ignore;
 		int use[3];
-	} COMBO3[4] = {
+	} COMBO3[] = {
 		{ 0, { 1, 2, 3 } },
 		{ 1, { 0, 2, 3 } },
 		{ 2, { 0, 1, 3 } },
 		{ 3, { 0, 1, 2 } }
 	};
 
-	buffermax = diskmax + 4 + 4 + 1 + 1;
+	buffermax = diskmax + LEV_MAX * 2 + 1 + 1;
 
 	buffer = malloc_nofail_vector_align(buffermax, block_size, &buffer_alloc);
 
@@ -440,7 +111,7 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 	}
 
 	/* compute the parity */
-	raid_gen(4, buffer, diskmax, block_size);
+	raid_gen(LEV_MAX, buffer, diskmax, block_size);
 
 	/* check 4data */
 	for(i=0;i<diskmax;++i) {
@@ -473,7 +144,7 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 						|| memcmp(buffer_test2, c, block_size) != 0
 						|| memcmp(buffer_test3, d, block_size) != 0
 					) {
-						fprintf(stderr, "Failed RAIDQP test\n");
+						fprintf(stderr, "Failed RAID test\n");
 						exit(EXIT_FAILURE);
 					}
 
@@ -491,12 +162,16 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 	for(i=0;i<diskmax;++i) {
 		for(j=i+1;j<diskmax;++j) {
 			for(k=j+1;k<diskmax;++k) {
-				for(n=0;n<4;++n) {
+				for(n=0;n<sizeof(COMBO3)/sizeof(COMBO3[0]);++n) {
 					unsigned char* a;
 					unsigned char* b;
 					unsigned char* c;
 					unsigned char* p;
 
+					/* this is the point where quad parity fails over RAID_QUAD_DATA_LIMIT data disks */
+					if (diskmax > RAID_QUAD_DATA_LIMIT && COMBO3[n].use[2] >= 3)
+						continue;
+
 					/* save */
 					a = buffer[i];
 					b = buffer[j];
@@ -520,7 +195,7 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 						|| memcmp(buffer_test1, b, block_size) != 0
 						|| memcmp(buffer_test2, c, block_size) != 0
 					) {
-						fprintf(stderr, "Failed RAIDQP test 1\n");
+						fprintf(stderr, "Failed RAID test\n");
 						exit(EXIT_FAILURE);
 					}
 
@@ -539,7 +214,7 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 	/* check 2data */
 	for(i=0;i<diskmax;++i) {
 		for(j=i+1;j<diskmax;++j) {
-			for(n=0;n<6;++n) {
+			for(n=0;n<sizeof(COMBO2)/sizeof(COMBO2[0]);++n) {
 				unsigned char* a;
 				unsigned char* b;
 				unsigned char* p;
@@ -568,7 +243,7 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 				if (memcmp(buffer_test0, a, block_size) != 0
 					|| memcmp(buffer_test1, b, block_size) != 0
 				) {
-					fprintf(stderr, "Failed RAIDQP test 1\n");
+					fprintf(stderr, "Failed RAID test\n");
 					exit(EXIT_FAILURE);
 				}
 
@@ -584,7 +259,7 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 	}
 
 	/* check 1data */
-	for(n=0;n<4;++n) {
+	for(n=0;n<sizeof(COMBO1)/sizeof(COMBO1[0]);++n) {
 		for(i=0;i<diskmax;++i) {
 			unsigned char* a;
 			unsigned char* p;
@@ -613,7 +288,7 @@ static void raidQPtest(unsigned diskmax, unsigned block_size)
 
 			/* check */
 			if (memcmp(buffer_test0, a, block_size) != 0) {
-				fprintf(stderr, "Failed RAIDQP test\n");
+				fprintf(stderr, "Failed RAID test\n");
 				exit(EXIT_FAILURE);
 			}
 
@@ -746,9 +421,6 @@ void selftest(int gui)
 
 	hashtest();
 	crc32ctest();
-	raid5test(33, 1024);
-	raid6test(33, 1024);
-	raidTPtest(33, 1024);
-	raidQPtest(21, 1024);
+	raidtest(24, 2048);
 }
 

From 81d3b2ad78ae60c9c85de1909701565b2b15e74d Mon Sep 17 00:00:00 2001
From: Andrea Mazzoleni <amadvance@gmail.com>
Date: Wed, 18 Sep 2013 23:16:39 +0200
Subject: [PATCH 4/4] Enforces the 21 data disk limit.

---
 raid.h     |  9 +++++++--
 selftest.c |  4 ++--
 state.c    | 12 +++++++++---
 3 files changed, 18 insertions(+), 7 deletions(-)

diff --git a/raid.h b/raid.h
index 14fb1194..9a7c74ad 100644
--- a/raid.h
+++ b/raid.h
@@ -22,9 +22,14 @@
 /* raid */
 
 /**
- * Max number of data disks for quad parity.
+ * Max number of data disks for RAID6.
  */
-#define RAID_QUAD_DATA_LIMIT 21
+#define RAID6_DATA_LIMIT 21
+
+/**
+ * Max number of data disks for RAIDQP.
+ */
+#define RAIDQP_DATA_LIMIT 21
 
 /**
  * Syndrome computation.
diff --git a/selftest.c b/selftest.c
index 3e2204b3..bbc99c66 100644
--- a/selftest.c
+++ b/selftest.c
@@ -168,8 +168,8 @@ static void raidtest(unsigned diskmax, unsigned block_size)
 					unsigned char* c;
 					unsigned char* p;
 
-					/* this is the point where quad parity fails over RAID_QUAD_DATA_LIMIT data disks */
-					if (diskmax > RAID_QUAD_DATA_LIMIT && COMBO3[n].use[2] >= 3)
+					/* this is the point where quad parity fails over RAIDQP_DATA_LIMIT data disks */
+					if (diskmax > RAIDQP_DATA_LIMIT && COMBO3[n].use[2] >= 3)
 						continue;
 
 					/* save */
diff --git a/state.c b/state.c
index 320df9e1..d6942f50 100644
--- a/state.c
+++ b/state.c
@@ -23,6 +23,7 @@
 #include "util.h"
 #include "parity.h"
 #include "cpu.h"
+#include "raid.h"
 
 const char* lev_name(unsigned l)
 {
@@ -176,9 +177,14 @@ static void state_config_check(struct snapraid_state* state, const char* path)
 			++diskcount;
 		}
 
-		if (diskcount > 255) {
-			/* RAID6 P/Q parity works for up to 255 drives, no more */
-			fprintf(stderr, "Too many disks. No more than 255.\n");
+		if (state->level >= 2 && diskcount > RAID6_DATA_LIMIT) {
+			/* RAID6 parity works for up to 255 drives, no more */
+			fprintf(stderr, "Too many disks for RAID6. No more than %u.\n", RAID6_DATA_LIMIT);
+			exit(EXIT_FAILURE);
+		}
+		if (state->level >= 4 && diskcount > RAIDQP_DATA_LIMIT) {
+			/* RAIDQP parity works for up to 21 drives, no more */
+			fprintf(stderr, "Too many disks for RAIDQP. No more than %u.\n", RAIDQP_DATA_LIMIT);
 			exit(EXIT_FAILURE);
 		}