diff --git a/check.c b/check.c
index c94e87dc..1ec85967 100644
--- a/check.c
+++ b/check.c
@@ -151,7 +151,7 @@ static int repair_step(struct snapraid_state* state, int rehash, unsigned pos, u
 		/* copy the redundancy to use */
 		memcpy(buffer[diskmax+1], buffer_qarity, state->block_size);
 
-		/* recover both data and p */
+		/* recover data */
 		raid6_recov_datap(buffer, diskmax, state->block_size, failed[failed_map[0]].index, buffer_zero);
 
 		hash_checked = 0; /* keep track if we check at least one block */
@@ -181,9 +181,12 @@ static int repair_step(struct snapraid_state* state, int rehash, unsigned pos, u
 		/* copy the redundancy to use */
 		memcpy(buffer[diskmax+1], buffer_qarity, state->block_size);
 
-		/* recover both data and p */
+		/* recover data */
 		raid6_recov_datap(buffer, diskmax, state->block_size, failed[failed_map[0]].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 (memcmp(buffer[diskmax], buffer_parity, state->block_size) == 0) {
 			/* recompute all the redundancy information */
diff --git a/raid.c b/raid.c
index bee35ea5..1dc349ce 100644
--- a/raid.c
+++ b/raid.c
@@ -33,12 +33,11 @@
 #include "cpu.h"
 #include "tables.h"
 
-/* internal forwarder */
-static void (*raid5_gen)(unsigned char** buffer, unsigned diskmax, unsigned size);
-static void (*raid6_gen)(unsigned char** buffer, unsigned diskmax, unsigned size);
+/****************************************************************************/
+/* internal computation */
 
 /*
- * 2-way unrolled portable integer math RAID-5 instruction set
+ * RAID5 unrolled-by-2 C implementation
  */
 void raid5_int32r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 {
@@ -49,8 +48,8 @@ void raid5_int32r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 	uint32_t wp0;
 	uint32_t wp1;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
 
 	for(d=0;d<size;d+=8) {
 		wp0 = *(uint32_t *)&buffer[z0][d+0*4];
@@ -65,9 +64,8 @@ void raid5_int32r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 }
 
 #if defined(__i386__) || defined(__x86_64__)
-
 /*
- * Unrolled-by-2 MMX implementation
+ * RAID5 unrolled-by-2 MMX implementation
  */
 void raid5_mmxr2(unsigned char** buffer, unsigned diskmax, unsigned size)
 {
@@ -75,8 +73,8 @@ void raid5_mmxr2(unsigned char** buffer, unsigned diskmax, unsigned size)
 	int z, z0;
 	unsigned d;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
 
 	for(d=0;d<size;d+=16) {
 		asm volatile("movq %0,%%mm0" : : "m" (buffer[z0][d]));
@@ -95,7 +93,7 @@ void raid5_mmxr2(unsigned char** buffer, unsigned diskmax, unsigned size)
 }
 
 /*
- * Unrolled-by-4 MMX implementation
+ * RAID5 unrolled-by-4 MMX implementation
  */
 void raid5_mmxr4(unsigned char** buffer, unsigned diskmax, unsigned size)
 {
@@ -103,8 +101,8 @@ void raid5_mmxr4(unsigned char** buffer, unsigned diskmax, unsigned size)
 	int z, z0;
 	unsigned d;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
 
 	for(d=0;d<size;d+=32) {
 		asm volatile("movq %0,%%mm0" : : "m" (buffer[z0][d]));
@@ -131,7 +129,7 @@ void raid5_mmxr4(unsigned char** buffer, unsigned diskmax, unsigned size)
 }
 
 /*
- * Unrolled-by-2 SSE2 implementation
+ * RAID5 unrolled-by-2 SSE2 implementation
  */
 void raid5_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 {
@@ -139,8 +137,8 @@ void raid5_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 	int z, z0;
 	unsigned d;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
 
 	for(d=0;d<size;d+=32) {
 		asm volatile("movdqa %0,%%xmm0" : : "m" (buffer[z0][d]));
@@ -159,7 +157,7 @@ void raid5_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 }
 
 /*
- * Unrolled-by-4 SSE2 implementation
+ * RAID5 unrolled-by-4 SSE2 implementation
  */
 void raid5_sse2r4(unsigned char** buffer, unsigned diskmax, unsigned size)
 {
@@ -167,8 +165,8 @@ void raid5_sse2r4(unsigned char** buffer, unsigned diskmax, unsigned size)
 	int z, z0;
 	unsigned d;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
 
 	for(d=0;d<size;d+=64) {
 		asm volatile("movdqa %0,%%xmm0" : : "m" (buffer[z0][d]));
@@ -197,7 +195,7 @@ void raid5_sse2r4(unsigned char** buffer, unsigned diskmax, unsigned size)
 
 #if defined(__x86_64__)
 /*
- * Unrolled-by-8 SSE2 implementation
+ * RAID5 unrolled-by-8 SSE2 implementation
  * Note that it uses all the 16 registers, meaning that x64 is required.
  */
 void raid5_sse2r8(unsigned char** buffer, unsigned diskmax, unsigned size)
@@ -206,8 +204,8 @@ void raid5_sse2r8(unsigned char** buffer, unsigned diskmax, unsigned size)
 	int z, z0;
 	unsigned d;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
 
 	for(d=0;d<size;d+=128) {
 		asm volatile("movdqa %0,%%xmm0" : : "m" (buffer[z0][d]));
@@ -286,7 +284,7 @@ static inline uint32_t MASK(uint32_t v)
 }
 
 /*
- * 2-way unrolled portable integer math RAID-6 instruction set
+ * RAID6 unrolled-by-2 C implementation
  */
 void raid6_int32r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 {
@@ -298,9 +296,9 @@ void raid6_int32r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 	uint32_t wd0, wq0, wp0, w10, w20;
 	uint32_t wd1, wq1, wp1, w11, w21;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
-	q = buffer[z0+2]; /* RS qarity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
 
 	for(d=0;d<size;d+=8) {
 		wq0 = wp0 = *(uint32_t *)&buffer[z0][d+0*4];
@@ -329,7 +327,6 @@ void raid6_int32r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 }
 
 #if defined(__i386__) || defined(__x86_64__)
-
 static const struct raid6_mmx_constants {
 	uint64_t x1d;
 } raid6_mmx_constants = {
@@ -337,7 +334,7 @@ static const struct raid6_mmx_constants {
 };
 
 /*
- * Unrolled-by-2 MMX implementation
+ * RAID6 unrolled-by-2 MMX implementation
  */
 void raid6_mmxr2(unsigned char** buffer, unsigned diskmax, unsigned size)
 {
@@ -346,13 +343,13 @@ void raid6_mmxr2(unsigned char** buffer, unsigned diskmax, unsigned size)
 	int z, z0;
 	unsigned d;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
-	q = buffer[z0+2]; /* RS qarity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
 
 	asm volatile("movq %0,%%mm0" : : "m" (raid6_mmx_constants.x1d));
-	asm volatile("pxor %mm5,%mm5"); /* Zero temp */
-	asm volatile("pxor %mm7,%mm7"); /* Zero temp */
+	asm volatile("pxor %mm5,%mm5");
+	asm volatile("pxor %mm7,%mm7");
 
 	for(d=0;d<size;d+=16) {
 		asm volatile("movq %0,%%mm2" : : "m" (buffer[z0][d]));
@@ -393,7 +390,7 @@ static const struct raid6_sse_constants {
 };
 
 /*
- * Unrolled-by-2 SSE2 implementation
+ * RAID6 unrolled-by-2 SSE2 implementation
  */
 void raid6_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 {
@@ -402,13 +399,13 @@ void raid6_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 	int z, z0;
 	unsigned d;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
-	q = buffer[z0+2]; /* RS qarity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
 
 	asm volatile("movdqa %0,%%xmm7" : : "m" (raid6_sse_constants.x1d[0]));
-	asm volatile("pxor %xmm4,%xmm4"); /* Zero temp */
-	asm volatile("pxor %xmm5,%xmm5"); /* Zero temp */
+	asm volatile("pxor %xmm4,%xmm4");
+	asm volatile("pxor %xmm5,%xmm5");
 
 	for(d=0;d<size;d+=32) {
 		asm volatile("movdqa %0,%%xmm0" : : "m" (buffer[z0][d])); /* P[0] */
@@ -445,7 +442,7 @@ void raid6_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
 
 #if defined(__x86_64__)
 /*
- * Unrolled-by-4 SSE2 implementation
+ * RAID6 unrolled-by-4 SSE2 implementation
  * Note that it uses all the 16 registers, meaning that x64 is required.
  */
 void raid6_sse2r4(unsigned char** buffer, unsigned diskmax, unsigned size)
@@ -455,15 +452,15 @@ void raid6_sse2r4(unsigned char** buffer, unsigned diskmax, unsigned size)
 	int z, z0;
 	unsigned d;
 
-	z0 = diskmax - 1; /* Highest data disk */
-	p = buffer[z0+1]; /* XOR parity */
-	q = buffer[z0+2]; /* RS qarity */
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
 
 	asm volatile("movdqa %0,%%xmm15" : : "m" (raid6_sse_constants.x1d[0]));
-	asm volatile("pxor %xmm8,%xmm8"); /* Zero temp */
-	asm volatile("pxor %xmm9,%xmm9"); /* Zero temp */
-	asm volatile("pxor %xmm10,%xmm10"); /* Zero temp */
-	asm volatile("pxor %xmm11,%xmm11"); /* Zero temp */
+	asm volatile("pxor %xmm8,%xmm8");
+	asm volatile("pxor %xmm9,%xmm9");
+	asm volatile("pxor %xmm10,%xmm10");
+	asm volatile("pxor %xmm11,%xmm11");
 
 	for(d=0;d<size;d+=64) {
 		asm volatile("movdqa %0,%%xmm0" : : "m" (buffer[z0][d])); /* P[0] */
@@ -522,36 +519,394 @@ void raid6_sse2r4(unsigned char** buffer, unsigned diskmax, unsigned size)
 }
 #endif
 
+/*
+ * RAIDTP unrolled-by-2 C implementation
+ */
+void raidTP_int32r2(unsigned char** buffer, unsigned diskmax, unsigned size)
+{
+	unsigned char* p;
+	unsigned char* q;
+	unsigned char* r;
+	int z, z0;
+	unsigned d;
+
+	uint32_t wd0, wr0, wq0, wp0, w10, w20;
+	uint32_t wd1, wr1, wq1, wp1, w11, w21;
+
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
+	r = buffer[diskmax+2];
+
+	for(d=0;d<size;d+=8) {
+		wr0 = wq0 = wp0 = *(uint32_t *)&buffer[z0][d+0*4];
+		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;
+		}
+		*(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;
+	}
+}
+
+#if defined(__i386__) || defined(__x86_64__)
+/*
+ * RAIDTP unrolled-by-1 MMX implementation
+ */
+void raidTP_mmxr1(unsigned char** buffer, unsigned diskmax, unsigned size)
+{
+	unsigned char* p;
+	unsigned char* q;
+	unsigned char* r;
+	int z, z0;
+	unsigned d;
+
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
+	r = buffer[diskmax+2];
+
+	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"); /* Q[0] */
+		for(z=z0-1;z>=0;--z) {
+			asm volatile("pcmpgtb %mm1,%mm4");
+			asm volatile("pcmpgtb %mm2,%mm5");
+			asm volatile("paddb %mm1,%mm1");
+			asm volatile("paddb %mm2,%mm2");
+			asm volatile("pand %mm7,%mm4");
+			asm volatile("pand %mm7,%mm5");
+			asm volatile("pxor %mm4,%mm1");
+			asm volatile("pxor %mm5,%mm2");
+			asm volatile("pcmpgtb %mm2,%mm6");
+			asm volatile("paddb %mm2,%mm2");
+			asm volatile("pand %mm7,%mm6");
+			asm volatile("pxor %mm6,%mm2");
+			asm volatile("movq %0,%%mm3" : : "m" (buffer[z][d]));
+			asm volatile("pxor %mm3,%mm0");
+			asm volatile("pxor %mm3,%mm1");
+			asm volatile("pxor %mm3,%mm2");
+			asm volatile("pxor %mm4,%mm4");
+			asm volatile("pxor %mm5,%mm5");
+			asm volatile("pxor %mm6,%mm6");
+		}
+		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("emms" : : : "memory");
+}
+#endif
+
+#if defined(__i386__) || defined(__x86_64__)
+/*
+ * RAIDTP unrolled-by-1 SSE2 implementation
+ */
+void raidTP_sse2r1(unsigned char** buffer, unsigned diskmax, unsigned size)
+{
+	unsigned char* p;
+	unsigned char* q;
+	unsigned char* r;
+	int z, z0;
+	unsigned d;
+
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
+	r = buffer[diskmax+2];
+
+	asm volatile("movdqa %0,%%xmm7" : : "m" (raid6_sse_constants.x1d[0]));
+	asm volatile("pxor %xmm4,%xmm4");
+	asm volatile("pxor %xmm5,%xmm5");
+	asm volatile("pxor %xmm6,%xmm6");
+
+	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] */
+		for(z=z0-1;z>=0;--z) {
+			asm volatile("pcmpgtb %xmm1,%xmm4");
+			asm volatile("pcmpgtb %xmm2,%xmm5");
+			asm volatile("paddb %xmm1,%xmm1");
+			asm volatile("paddb %xmm2,%xmm2");
+			asm volatile("pand %xmm7,%xmm4");
+			asm volatile("pand %xmm7,%xmm5");
+			asm volatile("pxor %xmm4,%xmm1");
+			asm volatile("pxor %xmm5,%xmm2");
+			asm volatile("pcmpgtb %xmm2,%xmm6");
+			asm volatile("paddb %xmm2,%xmm2");
+			asm volatile("pand %xmm7,%xmm6");
+			asm volatile("pxor %xmm6,%xmm2");
+			asm volatile("movdqa %0,%%xmm3" : : "m" (buffer[z][d]));
+			asm volatile("pxor %xmm3,%xmm0");
+			asm volatile("pxor %xmm3,%xmm1");
+			asm volatile("pxor %xmm3,%xmm2");
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pxor %xmm5,%xmm5");
+			asm volatile("pxor %xmm6,%xmm6");
+		}
+		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("sfence" : : : "memory");
+}
+#endif
+
+
+#if defined(__x86_64__)
+/*
+ * RAIDTP unrolled-by-2 SSE2 implementation
+ * Note that it uses all the 16 registers, meaning that x64 is required.
+ */
+void raidTP_sse2r2(unsigned char** buffer, unsigned diskmax, unsigned size)
+{
+	unsigned char* p;
+	unsigned char* q;
+	unsigned char* r;
+	int z, z0;
+	unsigned d;
+
+	z0 = diskmax - 1;
+	p = buffer[diskmax];
+	q = buffer[diskmax+1];
+	r = buffer[diskmax+2];
+
+	asm volatile("movdqa %0,%%xmm7" : : "m" (raid6_sse_constants.x1d[0]));
+	asm volatile("pxor %xmm4,%xmm4");
+	asm volatile("pxor %xmm5,%xmm5");
+	asm volatile("pxor %xmm6,%xmm6");
+	asm volatile("pxor %xmm12,%xmm12");
+	asm volatile("pxor %xmm13,%xmm13");
+	asm volatile("pxor %xmm14,%xmm14");
+
+	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] */
+		for(z=z0-1;z>=0;--z) {
+			asm volatile("pcmpgtb %xmm1,%xmm4");
+			asm volatile("pcmpgtb %xmm2,%xmm5");
+			asm volatile("pcmpgtb %xmm9,%xmm12");
+			asm volatile("pcmpgtb %xmm10,%xmm13");
+			asm volatile("paddb %xmm1,%xmm1");
+			asm volatile("paddb %xmm2,%xmm2");
+			asm volatile("paddb %xmm9,%xmm9");
+			asm volatile("paddb %xmm10,%xmm10");
+			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,%xmm2");
+			asm volatile("pxor %xmm12,%xmm9");
+			asm volatile("pxor %xmm13,%xmm10");
+			asm volatile("pcmpgtb %xmm2,%xmm6");
+			asm volatile("pcmpgtb %xmm10,%xmm14");
+			asm volatile("paddb %xmm2,%xmm2");
+			asm volatile("paddb %xmm10,%xmm10");
+			asm volatile("pand %xmm7,%xmm6");
+			asm volatile("pand %xmm7,%xmm14");
+			asm volatile("pxor %xmm6,%xmm2");
+			asm volatile("pxor %xmm14,%xmm10");
+			asm volatile("movdqa %0,%%xmm3" : : "m" (buffer[z][d]));
+			asm volatile("movdqa %0,%%xmm11" : : "m" (buffer[z][d+16]));
+			asm volatile("pxor %xmm3,%xmm0");
+			asm volatile("pxor %xmm3,%xmm1");
+			asm volatile("pxor %xmm3,%xmm2");
+			asm volatile("pxor %xmm11,%xmm8");
+			asm volatile("pxor %xmm11,%xmm9");
+			asm volatile("pxor %xmm11,%xmm10");
+			asm volatile("pxor %xmm4,%xmm4");
+			asm volatile("pxor %xmm5,%xmm5");
+			asm volatile("pxor %xmm6,%xmm6");
+			asm volatile("pxor %xmm12,%xmm12");
+			asm volatile("pxor %xmm13,%xmm13");
+			asm volatile("pxor %xmm14,%xmm14");
+		}
+		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("sfence" : : : "memory");
+}
+#endif
+
+/****************************************************************************/
+/* generic computation */
+
+/* internal forwarder */
+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);
+
+void raid_gen(unsigned level, unsigned char** buffer, unsigned diskmax, unsigned size)
+{
+	switch (level) {
+	case 1 : raid5_gen(buffer, diskmax, size); break;
+	case 2 : raid6_gen(buffer, diskmax, size); break;
+	case 3 : raidTP_gen(buffer, diskmax, size); break;
+	}
+}
+
+/****************************************************************************/
+/* recovering */
+
+/**
+ * GF a*b.
+ */
+static inline unsigned char mul(unsigned char a, unsigned char b)
+{
+	return gfmul[a][b];
+}
+
+/**
+ * GF 1/a.
+ * Not defined for a == 0.
+ */
+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 multiplication table.
+ */
+static inline const unsigned char* table(unsigned char a)
+{
+	return gfmul[a];
+}
+
 void raid5_recov_data(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila)
 {
-	unsigned char *p, *d;
+	unsigned char* p;
+	unsigned char* dp;
 
-	d = dptrs[faila];
+	dp = dptrs[faila];
 	p = dptrs[diskmax];
 
-	/* Compute syndrome using parity for the missing data page. */
+	/* compute syndrome using parity for the missing data page */
 	dptrs[faila] = p;
-	dptrs[diskmax] = d;
+	dptrs[diskmax] = dp;
 
 	raid5_gen(dptrs, diskmax, size);
 
-	/* Restore pointer table */
-	dptrs[faila] = d;
+	/* restore pointers */
+	dptrs[faila] = dp;
 	dptrs[diskmax] = p;
 }
 
 void raid6_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, unsigned char* zero)
 {
-	unsigned char *p, *q, *dp, *dq;
-	unsigned char px, qx, db;
-	const unsigned char *pbmul; /* P multiplier table for B data */
-	const unsigned char *qmul; /* Q multiplier table (for both) */
+	unsigned char* p;
+	unsigned char* dp;
+	unsigned char* q;
+	unsigned char* dq;
+	const unsigned char* pamul; /* P multipler to compute A */
+	const unsigned char* qamul; /* Q multipler to compute A */
 
 	p = dptrs[diskmax];
 	q = dptrs[diskmax+1];
 
-	/* Compute syndrome with zero for the missing data pages. */
-	/* Use the dead data pages as temporary storage for delta p and delta q. */
+	/* compute syndrome with zero for the missing data pages */
+	/* use the dead data pages as temporary storage for delta p and delta q. */
 	dp = dptrs[faila];
 	dptrs[faila] = zero;
 	dptrs[diskmax] = dp;
@@ -561,83 +916,367 @@ void raid6_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, i
 
 	raid6_gen(dptrs, diskmax, size);
 
-	/* Restore pointer table */
+	/* restore pointers */
 	dptrs[faila] = dp;
 	dptrs[failb] = dq;
 	dptrs[diskmax] = p;
 	dptrs[diskmax+1] = q;
 
-	/* Now, pick the proper data tables */
-	pbmul = raid6_gfmul[raid6_gfexi[failb-faila]];
-	qmul = raid6_gfmul[raid6_gfinv[raid6_gfexp[faila]^raid6_gfexp[failb]]];
+	/* select tables */
+	pamul = table( inv(pow2(failb-faila) ^ 1) );
+	qamul = table( inv(pow2(faila) ^ pow2(failb)) );
 
-	/* Now do it... */
 	while (size--) {
-		px = *p ^ *dp;
-		qx = qmul[*q ^ *dq];
-		*dq++ = db = pbmul[px] ^ qx; /* Reconstructed B */
-		*dp++ = db ^ px; /* Reconstructed A */
-		p++;
-		q++;
+		/* delta */
+		unsigned char pd = *p ^ *dp;
+		unsigned char qd = *q ^ *dq;
+
+		/* addends to reconstruct A */
+		unsigned char pm = pamul[pd];
+		unsigned char qm = qamul[qd];
+
+		/* reconstruct B */
+		unsigned char b = pm ^ qm;
+
+		/* reconstruct A */
+		unsigned char a = pd ^ b;
+
+		/* set */
+		*dp = a;
+		*dq = b;
+
+		++p;
+		++dp;
+		++q;
+		++dq;
 	}
 }
 
 void raid6_recov_datap(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, unsigned char* zero)
 {
-	unsigned char *p, *q, *dq;
-	const unsigned char* qmul; /* Q multiplier table */
+	unsigned char* q;
+	unsigned char* dq;
+	const unsigned char* qamul; /* Q multipler to compute A */
 
-	p = dptrs[diskmax];
 	q = dptrs[diskmax+1];
 
-	/* Compute syndrome with zero for the missing data page. */
-	/* Use the dead data page as temporary storage for delta q. */
+	/* compute syndrome with zero for the missing data page */
+	/* use the dead data page as temporary storage for delta q */
 	dq = dptrs[faila];
 	dptrs[faila] = zero;
 	dptrs[diskmax+1] = dq;
 
 	raid6_gen(dptrs, diskmax, size);
 
-	/* Restore pointer table */
+	/* restore pointers */
 	dptrs[faila] = dq;
 	dptrs[diskmax+1] = q;
 
-	/* Now, pick the proper data tables */
-	qmul = raid6_gfmul[raid6_gfinv[raid6_gfexp[faila]]];
+	/* select tables */
+	qamul = table( inv(pow2(faila)) );
 
-	/* Now do it... */
 	while (size--) {
-		*p++ ^= *dq = qmul[*q ^ *dq];
-		q++;
-		dq++;
+		/* delta */
+		unsigned char qd = *q ^ *dq;
+
+		/* addends to reconstruct A */
+		unsigned char qm = qamul[qd];
+
+		/* reconstruct A */
+		unsigned char a = qm;
+
+		/* set */
+		*dq = a;
+
+		++q;
+		++dq;
 	}
 }
 
-void raid_gen(unsigned level, unsigned char** buffer, unsigned diskmax, unsigned size)
+void raidTP_recov_datapq(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, unsigned char* zero)
 {
-	if (level == 1) {
-		raid5_gen(buffer, diskmax, size);
-	} else {
-		raid6_gen(buffer, diskmax, size);
+	unsigned char* r;
+	unsigned char* dr;
+	const unsigned char* ramul; /* R multiplier to compute A */
+
+	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[faila];
+	dptrs[faila] = zero;
+	dptrs[diskmax+2] = dr;
+
+	raidTP_gen(dptrs, diskmax, size);
+
+	/* restore pointers */
+	dptrs[faila] = dr;
+	dptrs[diskmax+2] = r;
+
+	/* select tables */
+	ramul = table( inv(pow4(faila)) );
+
+	while (size--) {
+		/* delta */
+		unsigned char rd = *r ^ *dr;
+
+		/* addends to reconstruct A */
+		unsigned char rm = ramul[rd];
+
+		/* reconstruct A */
+		unsigned char a = rm;
+
+		/* set */
+		*dr = a;
+
+		++r;
+		++dr;
 	}
 }
 
+void raidTP_recov_2dataq(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, unsigned char* zero)
+{
+	unsigned char* p;
+	unsigned char* dp;
+	unsigned char* r;
+	unsigned char* dr;
+	const unsigned char* pbmul; /* Q multipler to compute B */
+	const unsigned char* rbmul; /* R multipler to compute B */
+	unsigned char denominator;
+
+	p = dptrs[diskmax];
+	r = dptrs[diskmax+2];
+
+	/* compute syndrome with zero for the missing data page */
+	/* use the dead data page as temporary storage for delta p and r */
+	dp = dptrs[faila];
+	dr = dptrs[failb];
+	dptrs[faila] = zero;
+	dptrs[failb] = zero;
+	dptrs[diskmax] = dp;
+	dptrs[diskmax+2] = dr;
+
+	raidTP_gen(dptrs, diskmax, size);
+
+	/* restore pointers */
+	dptrs[faila] = dp;
+	dptrs[failb] = dr;
+	dptrs[diskmax] = p;
+	dptrs[diskmax+2] = r;
+
+	/* select tables */
+	denominator = pow4(failb-faila) ^ 1;
+	pbmul = table( inv(denominator) );
+	rbmul = table( inv(mul(denominator, pow4(faila))) );
+
+	while (size--) {
+		/* delta */
+		unsigned char pd = *p ^ *dp;
+		unsigned char rd = *r ^ *dr;
+
+		/* addends to reconstruct B */
+		unsigned char pbm = pbmul[pd];
+		unsigned char rbm = rbmul[rd];
+
+		/* reconstruct B */
+		unsigned char b = pbm ^ rbm;
+
+		/* reconstruct A */
+		unsigned char a = pd ^ b;
+
+		/* set */
+		*dp = a;
+		*dr = b;
+
+		++p;
+		++dp;
+		++r;
+		++dr;
+	}
+}
+
+void raidTP_recov_2datap(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, unsigned char* zero)
+{
+	unsigned char* q;
+	unsigned char* dq;
+	unsigned char* r;
+	unsigned char* dr;
+	const unsigned char* qbmul; /* Q multipler to compute B */
+	const unsigned char* rbmul; /* R multipler to compute B */
+	const unsigned char* qamul; /* Q multipler to compute A */
+	const unsigned char* bamul; /* B multipler to compute A */
+	unsigned char denominator;
+
+	q = dptrs[diskmax+1];
+	r = dptrs[diskmax+2];
+
+	/* compute syndrome with zero for the missing data page */
+	/* use the dead data page as temporary storage for delta q and r */
+	dq = dptrs[faila];
+	dr = dptrs[failb];
+	dptrs[faila] = zero;
+	dptrs[failb] = zero;
+	dptrs[diskmax+1] = dq;
+	dptrs[diskmax+2] = dr;
+
+	raidTP_gen(dptrs, diskmax, size);
+
+	/* restore pointers */
+	dptrs[faila] = dq;
+	dptrs[failb] = dr;
+	dptrs[diskmax+1] = q;
+	dptrs[diskmax+2] = r;
+
+	/* select tables */
+	denominator = pow2(failb-faila) ^ pow4(failb-faila);
+	qbmul = table( inv(mul(denominator, pow2(faila))) );
+	rbmul = table( inv(mul(denominator, pow4(faila))) );
+	qamul = table( inv(pow2(faila)) );
+	bamul = table( pow2(failb-faila) );
+
+	while (size--) {
+		/* delta */
+		unsigned char qd = *q ^ *dq;
+		unsigned char rd = *r ^ *dr;
+
+		/* addends to reconstruct B */
+		unsigned char qbm = qbmul[qd];
+		unsigned char rbm = rbmul[rd];
+
+		/* reconstruct B */
+		unsigned char b = qbm ^ rbm;
+
+		/* addends to reconstruct A */
+		unsigned char qam = qamul[qd];
+		unsigned char bam = bamul[b];
+
+		/* reconstruct A */
+		unsigned char a = qam ^ bam;
+
+		/* set */
+		*dq = a;
+		*dr = b;
+
+		++q;
+		++dq;
+		++r;
+		++dr;
+	}
+}
+
+void raidTP_recov_3data(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, int failc, unsigned char* zero)
+{
+	unsigned char* p;
+	unsigned char* dp;
+	unsigned char* q;
+	unsigned char* dq;
+	unsigned char* r;
+	unsigned char* dr;
+	const unsigned char* pcmul; /* P multipler to compute C */
+	const unsigned char* qcmul; /* Q multipler to compute C */
+	const unsigned char* rcmul; /* R multipler to compute C */
+	const unsigned char* pbmul; /* P multipler to compute B */
+	const unsigned char* qbmul; /* Q multipler to compute B */
+	const unsigned char* cbmul; /* C multipler to compute B */
+	unsigned char denominator;
+
+	p = dptrs[diskmax];
+	q = dptrs[diskmax+1];
+	r = dptrs[diskmax+2];
+
+	/* 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[faila];
+	dq = dptrs[failb];
+	dr = dptrs[failc];
+	dptrs[faila] = zero;
+	dptrs[failb] = zero;
+	dptrs[failc] = zero;
+	dptrs[diskmax] = dp;
+	dptrs[diskmax+1] = dq;
+	dptrs[diskmax+2] = dr;
+
+	raidTP_gen(dptrs, diskmax, size);
+
+	/* restore pointers */
+	dptrs[faila] = dp;
+	dptrs[failb] = dq;
+	dptrs[failc] = dr;
+	dptrs[diskmax] = p;
+	dptrs[diskmax+1] = q;
+	dptrs[diskmax+2] = r;
+
+	/* select tables */
+	denominator = mul(pow2(failc-faila) ^ 1, inv(pow2(failb-faila) ^ 1)) ^ mul(pow4(failc-faila) ^ 1, inv(pow4(failb-faila) ^ 1));
+	pcmul = table( inv(mul(denominator, pow2(failb-faila) ^ 1)) ^ inv(mul(denominator, pow4(failb-faila) ^ 1)) );
+	qcmul = table( inv(mul(mul(denominator, pow2(failb-faila) ^ 1), pow2(faila))) );
+	rcmul = table( inv(mul(mul(denominator, pow4(failb-faila) ^ 1), pow4(faila))) );
+	pbmul = table( inv(pow2(failb-faila) ^ 1) );
+	qbmul = table( inv(mul(pow2(faila), pow2(failb-faila) ^ 1)) );
+	cbmul = table( mul(pow2(failc-faila) ^ 1, inv(pow2(failb-faila) ^ 1)) );
+
+	while (size--) {
+		/* delta */
+		unsigned char pd = *p ^ *dp;
+		unsigned char qd = *q ^ *dq;
+		unsigned char rd = *r ^ *dr;
+
+		/* addends to reconstruct C */
+		unsigned char pcm = pcmul[pd];
+		unsigned char qcm = qcmul[qd];
+		unsigned char rcm = rcmul[rd];
+
+		/* reconstruct C */
+		unsigned char c = pcm ^ qcm ^ rcm;
+
+		/* addends to reconstruct B */
+		unsigned char pbm = pbmul[pd];
+		unsigned char qbm = qbmul[qd];
+		unsigned char cbm = cbmul[c];
+
+		/* reconstruct B */
+		unsigned char b = pbm ^ qbm ^ cbm;
+
+		/* reconstruct A */
+		unsigned char a = pd ^ b ^ c;
+
+		/* set */
+		*dp = a;
+		*dq = b;
+		*dr = c;
+
+		++p;
+		++dp;
+		++q;
+		++dq;
+		++r;
+		++dr;
+	}
+}
+
+/****************************************************************************/
+/* init/done */
+
 void raid_init(void)
 {
 	raid5_gen = raid5_int32r2;
 	raid6_gen = raid6_int32r2;
+	raidTP_gen = raidTP_int32r2;
 #if defined(__i386__) || defined(__x86_64__)
 	if (cpu_has_mmx()) {
 		raid5_gen = raid5_mmxr4;
 		raid6_gen = raid6_mmxr2;
+		raidTP_gen = raidTP_mmxr1;
 	}
 	if (cpu_has_sse2()) {
 #if defined(__x86_64__)
 		raid5_gen = raid5_sse2r8;
 		raid6_gen = raid6_sse2r4;
+		raidTP_gen = raidTP_sse2r2;
 #else
 		raid5_gen = raid5_sse2r4;
 		raid6_gen = raid6_sse2r2;
+		raidTP_gen = raidTP_sse2r1;
 #endif
 	}
 #endif
diff --git a/raid.h b/raid.h
index 7df08a37..f3df6e4f 100644
--- a/raid.h
+++ b/raid.h
@@ -37,10 +37,34 @@ void raid5_recov_data(unsigned char** dptrs, unsigned diskmax, unsigned size, in
 void raid6_recov_2data(unsigned char** dptrs, unsigned diskmax, unsigned size, int faila, int failb, unsigned char* zero);
 
 /*
- * Recover failure of one data block plus the P block.
+ * 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);
+
 /**
  * Initializes the RAID system.
  */
@@ -59,6 +83,10 @@ void raid6_int32r2(unsigned char** dptr, unsigned diskmax, unsigned size);
 void raid6_mmxr2(unsigned char** dptr, unsigned diskmax, unsigned size);
 void raid6_sse2r2(unsigned char** dptr, unsigned diskmax, unsigned size);
 void raid6_sse2r4(unsigned char** buffer, unsigned diskmax, unsigned size);
+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);
 
 #endif
 
diff --git a/selftest.c b/selftest.c
index 84e982e4..412f0cb7 100644
--- a/selftest.c
+++ b/selftest.c
@@ -47,7 +47,7 @@ static void raid5test(unsigned diskmax, unsigned block_size)
 {
 	void* buffer_alloc;
 	unsigned char** buffer;
-	unsigned char* buffer_save;
+	unsigned char* buffer_test;
 	unsigned buffermax;
 	unsigned i, j;
 
@@ -55,7 +55,7 @@ static void raid5test(unsigned diskmax, unsigned block_size)
 
 	buffer = malloc_nofail_vector_align(buffermax, block_size, &buffer_alloc);
 
-	buffer_save = buffer[buffermax-1];
+	buffer_test = buffer[buffermax-1];
 
 	/* fill with random */
 	for(i=0;i<diskmax;++i) {
@@ -63,20 +63,26 @@ static void raid5test(unsigned diskmax, unsigned block_size)
 			buffer[i][j] = rand();
 	}
 
+	/* compute the parity */
 	raid_gen(1, buffer, diskmax, block_size);
 
 	for(i=0;i<diskmax;++i) {
-		/* save the correct one */
-		memcpy(buffer_save, buffer[i], block_size);
+		unsigned char* a;
 
-		/* destroy it */
-		memset(buffer[i], 0x55, block_size);
+		/* 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_save, buffer[i], block_size) != 0) {
+		if (memcmp(buffer_test, buffer[i], block_size) != 0) {
 			fprintf(stderr, "Failed RAID5 test\n");
 			exit(EXIT_FAILURE);
 		}
@@ -90,20 +96,23 @@ static void raid6test(unsigned diskmax, unsigned block_size)
 {
 	void* buffer_alloc;
 	unsigned char** buffer;
-	unsigned char* buffer_save0;
-	unsigned char* buffer_save1;
+	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;
+	buffermax = diskmax + 2 + 2 + 1 + 1;
 
 	buffer = malloc_nofail_vector_align(buffermax, block_size, &buffer_alloc);
 
-	buffer_save0 = buffer[buffermax-3];
-	buffer_save1 = buffer[buffermax-2];
-	buffer_zero = buffer[buffermax-1];
+	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) {
@@ -111,25 +120,33 @@ static void raid6test(unsigned diskmax, unsigned block_size)
 			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) {
-			/* save the correct ones */
-			memcpy(buffer_save0, buffer[i], block_size);
-			memcpy(buffer_save1, buffer[j], block_size);
+			unsigned char* a;
+			unsigned char* b;
 
-			/* destroy it */
-			memset(buffer[i], 0x55, block_size);
-			memset(buffer[j], 0x55, block_size);
+			/* 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_save0, buffer[i], block_size) != 0
-				|| memcmp(buffer_save1, buffer[j], block_size) != 0
+			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);
@@ -139,24 +156,29 @@ static void raid6test(unsigned diskmax, unsigned block_size)
 
 	/* check datap */
 	for(i=0;i<diskmax;++i) {
-		/* the other is the parity */
+		unsigned char* a;
+		unsigned char* p;
+
+		/* other is p */
 		j = diskmax;
 
-		/* save the correct ones */
-		memcpy(buffer_save0, buffer[i], block_size);
-		memcpy(buffer_save1, buffer[j], block_size);
+		/* save */
+		a = buffer[i];
+		p = buffer[j];
 
-		/* destroy it */
-		memset(buffer[i], 0x55, block_size);
-		memset(buffer[j], 0x55, block_size);
+		/* 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_save0, buffer[i], block_size) != 0
-			|| memcmp(buffer_save1, buffer[j], block_size) != 0
-		) {
+		if (memcmp(buffer_test0, buffer[i], block_size) != 0) {
 			fprintf(stderr, "Failed RAID6 test\n");
 			exit(EXIT_FAILURE);
 		}
@@ -166,6 +188,192 @@ static void raid6test(unsigned diskmax, unsigned block_size)
 	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);
+}
+
 #define HASH_TEST_MAX 512 /* tests are never longer than 512 bytes */
 
 static void hashtest(void)
@@ -201,7 +409,7 @@ static void hashtest(void)
 		memcpy(buffer_aligned, TEST_MURMUR3[i].data, TEST_MURMUR3[i].len);
 		memhash(HASH_MURMUR3, seed_aligned, digest, buffer_aligned, TEST_MURMUR3[i].len);
 		if (memcmp(digest, TEST_MURMUR3[i].digest, HASH_SIZE) != 0) {
-			fprintf(stderr, "Failed Murmur3 test vector\n");
+			fprintf(stderr, "Failed Murmur3 test\n");
 			exit(EXIT_FAILURE);
 		}
 	}
@@ -211,7 +419,7 @@ static void hashtest(void)
 		memcpy(buffer_aligned, TEST_SPOOKY2[i].data, TEST_SPOOKY2[i].len);
 		memhash(HASH_SPOOKY2, seed_aligned, digest, buffer_aligned, TEST_SPOOKY2[i].len);
 		if (memcmp(digest, TEST_SPOOKY2[i].digest, HASH_SIZE) != 0) {
-			fprintf(stderr, "Failed Spooky2 test vector %u\n", i);
+			fprintf(stderr, "Failed Spooky2 test\n");
 			exit(EXIT_FAILURE);
 		}
 	}
@@ -258,7 +466,7 @@ static void crc32ctest(void)
 		uint32_t digest;
 		digest = crc32c(0, (const unsigned char*)TEST_CRC32C[i].data, TEST_CRC32C[i].len);
 		if (digest != TEST_CRC32C[i].digest) {
-			fprintf(stderr, "Failed CRC32C test vector %u\n", i);
+			fprintf(stderr, "Failed CRC32C test\n");
 			exit(EXIT_FAILURE);
 		}
 	}
@@ -266,8 +474,6 @@ static void crc32ctest(void)
 
 void selftest(int gui)
 {
-	unsigned i;
-
 	if (gui) {
 		fprintf(stdlog, "selftest:\n");
 		fflush(stdlog);
@@ -283,10 +489,8 @@ void selftest(int gui)
 
 	hashtest();
 	crc32ctest();
-
-	for(i=1;i<=33;++i) {
-		raid5test(i, 2048);
-		raid6test(i, 2048);
-	}
+	raid5test(33, 2048);
+	raid6test(33, 2048);
+	raidTPtest(33, 2048);
 }
 
diff --git a/speed.c b/speed.c
index 1b6eab2a..c634e3a1 100644
--- a/speed.c
+++ b/speed.c
@@ -155,12 +155,24 @@ void speed(void)
 	else
 		printf("int32x2");
 	printf("\n");
+	printf("Expected fastest RAIDTP is ");
+	if (cpu_has_sse2())
+#if defined(__x86_64__)
+		printf("sse2x2");
+#else
+		printf("sse2x1");
+#endif
+	else if (cpu_has_mmx())
+		printf("mmxx1");
+	else
+		printf("int32x2");
+	printf("\n");
 #endif
 	printf("If these expectations are false, please report it in the SnapRAID forum\n");
 
 	printf("\n");
 
-	printf("Speed test with %d disk and %d buffer size, for a total of %u KiB...\n", diskmax, block_size, (diskmax + 2) * block_size / 1024);
+	printf("Speed test with %d disk and %d buffer size, for a total of %u KiB...\n", diskmax, block_size, (diskmax + 3) * block_size / 1024);
 
 	SPEED_START {
 		for(j=0;j<diskmax;++j) {
@@ -281,6 +293,38 @@ void speed(void)
 	}
 #endif
 
+	SPEED_START {
+		raidTP_int32r2(buffer, diskmax, block_size);
+	} SPEED_STOP
+
+	printf("RAIDTP int32x2 %"PRIu64" [MB/s]\n", ds / dt);
+
+#if defined(__i386__) || defined(__x86_64__)
+	if (cpu_has_mmx()) {
+		SPEED_START {
+			raidTP_mmxr1(buffer, diskmax, block_size);
+		} SPEED_STOP
+
+		printf("RAIDTP mmxr1 %"PRIu64" [MB/s]\n", ds / dt);
+	}
+
+	if (cpu_has_sse2()) {
+		SPEED_START {
+			raidTP_sse2r1(buffer, diskmax, block_size);
+		} SPEED_STOP
+
+		printf("RAIDTP sse2x1 %"PRIu64" [MB/s]\n", ds / dt);
+
+#if defined(__x86_64__)
+		SPEED_START {
+			raidTP_sse2r2(buffer, diskmax, block_size);
+		} SPEED_STOP
+
+		printf("RAIDTP sse2x2 %"PRIu64" [MB/s]\n", ds / dt);
+#endif
+	}
+#endif
+
 	free(buffer_alloc);
 	free(buffer);
 }
diff --git a/tables.h b/tables.h
index ad873ac3..205a6a85 100644
--- a/tables.h
+++ b/tables.h
@@ -15,7 +15,7 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-const unsigned char  __attribute__((aligned(256))) raid6_gfmul[256][256] =
+const unsigned char  __attribute__((aligned(256))) gfmul[256][256] =
 {
 	{
 		0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
@@ -8723,7 +8723,7 @@ const unsigned char  __attribute__((aligned(256))) raid6_gfmul[256][256] =
 	},
 };
 
-const unsigned char __attribute__((aligned(256))) raid6_gfexp[256] =
+const unsigned char __attribute__((aligned(256))) gfexp2[256] =
 {
 	0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80,
 	0x1d, 0x3a, 0x74, 0xe8, 0xcd, 0x87, 0x13, 0x26,
@@ -8757,10 +8757,49 @@ const unsigned char __attribute__((aligned(256))) raid6_gfexp[256] =
 	0xf7, 0xf3, 0xfb, 0xeb, 0xcb, 0x8b, 0x0b, 0x16,
 	0x2c, 0x58, 0xb0, 0x7d, 0xfa, 0xe9, 0xcf, 0x83,
 	0x1b, 0x36, 0x6c, 0xd8, 0xad, 0x47, 0x8e, 0x00,
+	/* note that the last element is not significative */
 };
 
-const unsigned char __attribute__((aligned(256))) raid6_gfinv[256] =
+const unsigned char __attribute__((aligned(256))) gfexp4[256] =
 {
+	0x01, 0x04, 0x10, 0x40, 0x1d, 0x74, 0xcd, 0x13,
+	0x4c, 0x2d, 0xb4, 0xea, 0x8f, 0x06, 0x18, 0x60,
+	0x9d, 0x4e, 0x25, 0x94, 0x6a, 0xb5, 0xee, 0x9f,
+	0x46, 0x05, 0x14, 0x50, 0x5d, 0x69, 0xb9, 0xde,
+	0x5f, 0x61, 0x99, 0x5e, 0x65, 0x89, 0x1e, 0x78,
+	0xfd, 0xd3, 0x6b, 0xb1, 0xfe, 0xdf, 0x5b, 0x71,
+	0xd9, 0x43, 0x11, 0x44, 0x0d, 0x34, 0xd0, 0x67,
+	0x81, 0x3e, 0xf8, 0xc7, 0x3b, 0xec, 0x97, 0x66,
+	0x85, 0x2e, 0xb8, 0xda, 0x4f, 0x21, 0x84, 0x2a,
+	0xa8, 0x9a, 0x52, 0x55, 0x49, 0x39, 0xe4, 0xb7,
+	0xe6, 0xbf, 0xc6, 0x3f, 0xfc, 0xd7, 0x7b, 0xf1,
+	0xe3, 0xab, 0x96, 0x62, 0x95, 0x6e, 0xa5, 0xae,
+	0x82, 0x32, 0xc8, 0x07, 0x1c, 0x70, 0xdd, 0x53,
+	0x51, 0x59, 0x79, 0xf9, 0xc3, 0x2b, 0xac, 0x8a,
+	0x12, 0x48, 0x3d, 0xf4, 0xf7, 0xfb, 0xcb, 0x0b,
+	0x2c, 0xb0, 0xfa, 0xcf, 0x1b, 0x6c, 0xad, 0x8e,
+	0x02, 0x08, 0x20, 0x80, 0x3a, 0xe8, 0x87, 0x26,
+	0x98, 0x5a, 0x75, 0xc9, 0x03, 0x0c, 0x30, 0xc0,
+	0x27, 0x9c, 0x4a, 0x35, 0xd4, 0x77, 0xc1, 0x23,
+	0x8c, 0x0a, 0x28, 0xa0, 0xba, 0xd2, 0x6f, 0xa1,
+	0xbe, 0xc2, 0x2f, 0xbc, 0xca, 0x0f, 0x3c, 0xf0,
+	0xe7, 0xbb, 0xd6, 0x7f, 0xe1, 0xa3, 0xb6, 0xe2,
+	0xaf, 0x86, 0x22, 0x88, 0x1a, 0x68, 0xbd, 0xce,
+	0x1f, 0x7c, 0xed, 0x93, 0x76, 0xc5, 0x33, 0xcc,
+	0x17, 0x5c, 0x6d, 0xa9, 0x9e, 0x42, 0x15, 0x54,
+	0x4d, 0x29, 0xa4, 0xaa, 0x92, 0x72, 0xd5, 0x73,
+	0xd1, 0x63, 0x91, 0x7e, 0xe5, 0xb3, 0xf6, 0xff,
+	0xdb, 0x4b, 0x31, 0xc4, 0x37, 0xdc, 0x57, 0x41,
+	0x19, 0x64, 0x8d, 0x0e, 0x38, 0xe0, 0xa7, 0xa6,
+	0xa2, 0xb2, 0xf2, 0xef, 0x9b, 0x56, 0x45, 0x09,
+	0x24, 0x90, 0x7a, 0xf5, 0xf3, 0xeb, 0x8b, 0x16,
+	0x58, 0x7d, 0xe9, 0x83, 0x36, 0xd8, 0x47, 0x00,
+	/* note that the last element is not significative */
+};
+
+const unsigned char __attribute__((aligned(256))) gfinv[256] =
+{
+	/* note that the first element is not significative */
 	0x00, 0x01, 0x8e, 0xf4, 0x47, 0xa7, 0x7a, 0xba,
 	0xad, 0x9d, 0xdd, 0x98, 0x3d, 0xaa, 0x5d, 0x96,
 	0xd8, 0x72, 0xc0, 0x58, 0xe0, 0x3e, 0x4c, 0x66,
@@ -8795,39 +8834,3 @@ const unsigned char __attribute__((aligned(256))) raid6_gfinv[256] =
 	0x42, 0xd4, 0xe8, 0x75, 0x7f, 0xff, 0x7e, 0xfd,
 };
 
-const unsigned char __attribute__((aligned(256))) raid6_gfexi[256] =
-{
-	0x00, 0xf4, 0xa7, 0x9d, 0x72, 0xed, 0x5f, 0x54,
-	0xa0, 0xa9, 0xfb, 0x74, 0x85, 0xe4, 0xc0, 0x8a,
-	0x67, 0xdc, 0x26, 0x63, 0xf2, 0xe9, 0xb4, 0xd2,
-	0x02, 0x8e, 0xba, 0xaa, 0xde, 0x89, 0x57, 0x59,
-	0xac, 0x2c, 0x93, 0x09, 0x60, 0x1e, 0xbc, 0xa4,
-	0xb3, 0x43, 0xeb, 0xda, 0x6a, 0x12, 0xc7, 0x39,
-	0x04, 0x69, 0x47, 0x98, 0x3e, 0xd0, 0x23, 0x82,
-	0x97, 0x7b, 0x52, 0xf6, 0xc5, 0x2f, 0x91, 0x1c,
-	0x37, 0xb6, 0x24, 0x35, 0x0b, 0x33, 0x41, 0x9a,
-	0xb9, 0x5a, 0x49, 0x5d, 0x2a, 0x0c, 0x77, 0xe7,
-	0x7f, 0xe2, 0xc8, 0x07, 0xef, 0xd6, 0x87, 0xfe,
-	0xfd, 0x14, 0x19, 0x9f, 0xcb, 0xaf, 0x28, 0xf0,
-	0x10, 0xbe, 0xf8, 0xb0, 0xd8, 0x51, 0x4e, 0x31,
-	0x3d, 0x80, 0xc3, 0x8d, 0x71, 0x95, 0x17, 0x3a,
-	0x1b, 0x4b, 0xe1, 0xcd, 0xd4, 0x6d, 0xa3, 0x44,
-	0xcf, 0x79, 0x21, 0x65, 0x0f, 0x6f, 0x4d, 0x7d,
-	0x7c, 0x4c, 0x6e, 0x0e, 0x64, 0x20, 0x78, 0xce,
-	0x45, 0xa2, 0x6c, 0xd5, 0xcc, 0xe0, 0x4a, 0x1a,
-	0x3b, 0x16, 0x94, 0x70, 0x8c, 0xc2, 0x81, 0x3c,
-	0x30, 0x4f, 0x50, 0xd9, 0xb1, 0xf9, 0xbf, 0x11,
-	0xf1, 0x29, 0xae, 0xca, 0x9e, 0x18, 0x15, 0xfc,
-	0xff, 0x86, 0xd7, 0xee, 0x06, 0xc9, 0xe3, 0x7e,
-	0xe6, 0x76, 0x0d, 0x2b, 0x5c, 0x48, 0x5b, 0xb8,
-	0x9b, 0x40, 0x32, 0x0a, 0x34, 0x25, 0xb7, 0x36,
-	0x1d, 0x90, 0x2e, 0xc4, 0xf7, 0x53, 0x7a, 0x96,
-	0x83, 0x22, 0xd1, 0x3f, 0x99, 0x46, 0x68, 0x05,
-	0x38, 0xc6, 0x13, 0x6b, 0xdb, 0xea, 0x42, 0xb2,
-	0xa5, 0xbd, 0x1f, 0x61, 0x08, 0x92, 0x2d, 0xad,
-	0x58, 0x56, 0x88, 0xdf, 0xab, 0xbb, 0x8f, 0x03,
-	0xd3, 0xb5, 0xe8, 0xf3, 0x62, 0x27, 0xdd, 0x66,
-	0x8b, 0xc1, 0xe5, 0x84, 0x75, 0xfa, 0xa8, 0xa1,
-	0x55, 0x5e, 0xec, 0x73, 0x9c, 0xa6, 0xf5, 0x01,
-};
-