// SPDX-License-Identifier: GPL-2.0-or-later // Copyright (C) 2013 Andrea Mazzoleni #include "internal.h" #include "gf.h" /* * Extended Cauchy Matrix Construction for MDS RAID in GF(2^8) * * The Cauchy RAID modes operate in GF(2^8) and support up to six parity * levels. RAID_MODE_CAUCHY_RAID uses polynomial 0x11d, while * RAID_MODE_CAUCHY_AES uses polynomial 0x11b. * * For RAID5 and RAID6, it follows the method described in H. Peter Anvin's * paper "The mathematics of RAID-6" [1]. Please refer to this paper for a * complete explanation. * * Triple parity was first evaluated using an extension of the same approach, * with additional parity coefficients set as powers of 2^-1, defined as: * * P = sum(Di) * Q = sum(2^i * Di) * R = sum(2^-i * Di) for 0 <= i < N * * This approach works efficiently for triple parity because multiplications * and divisions by 2 in GF(2^8) are very fast. It is similar to ZFS RAIDZ3, * which uses powers of 4 instead of 2^-1. * * Unfortunately, this method does not extend beyond triple parity because * no choice of power coefficients guarantees solvable equations for all * combinations of missing disks. This is expected: a Vandermonde matrix, * used in this method, does not ensure that all submatrices are nonsingular * [2, Chap. 11, Problem 7], which is required for an MDS (Maximum Distance * Separable) code [2, Chap. 11, Theorem 8]. * * To overcome this limitation, a Cauchy matrix [3][4] is used for parity * computation. A Cauchy matrix ensures that all square submatrices are * nonsingular, so the linear system is always solvable for any combination * of missing disks. This guarantees an MDS code. * * RAID_MODE_CAUCHY_RAID uses polynomial 0x11d and the existing primitive * generator g=2 construction described below. * * RAID_MODE_CAUCHY_AES instead uses polynomial 0x11b and a G23 Extended * Cauchy construction. In the AES field, element 2 has multiplicative order * 51 while element 3 is primitive. Therefore H=<2> has five cosets in * GF(256)*. The Q coefficients q[0..254] enumerate those cosets, starting * with q[0]=1. Within a coset q advances by multiplication by 2. The four * transitions between cosets advance by multiplication by 3, after the * coefficients indexed 50, 101, 152, and 203. * * The AES data X set consumes inverses of the G23 sequence from the beginning: * * X[i] = inv(q[i]) * * The Cauchy Y set starts with zero and consumes inverses backward from the * end of the sequence: * * Y[0] = 0 * Y[j] = inv(q[255-j]) for j > 0 * * Thus the Y values associated with R, S, T, and U are inv(q[254]), * inv(q[253]), inv(q[252]), and inv(q[251]), respectively. For PARITY=n and * DISK=257-n, increasing parity removes one X element from the tail and * appends it to Y, leaving the existing rows unchanged on the remaining * columns. X and Y are disjoint and together contain all 256 field elements, so * C[j,i]=1/(X[i]+Y[j]) defines valid Cauchy rows. Adding P as an all-ones row * forms an Extended Cauchy matrix. Q is the row for Y[0]=0 and therefore * Q=sum(q[i]*Di). Rows R through U are scaled so their first coefficient is 1. * This nonzero row scaling preserves the MDS property, and the resulting 6x251 * matrix is MDS. * * The numerical examples below refer to RAID_MODE_CAUCHY_RAID, * using polynomial 0x11d and generator g = 2. * * How the RAID/g=2 matrix is obtained: * * The matrix is constructed to match Linux RAID coefficients for the * first two rows, while ensuring that all square submatrices are nonsingular. * * 1. Start by forming a Cauchy matrix with elements defined as 1/(x_i + y_j), * where all x_i and y_j are distinct elements (the textbook definition * of a Cauchy matrix). * * 2. For the first row (j=0), set x_i = g^-i and y_0 = 0, resulting in: * * row j=0 -> 1/(x_i + y_0) = 1/(g^-i + 0) = g^i * * which reproduces the RAID-6 coefficients. * * 3. For subsequent rows (j>0), set y_j = g^j, yielding: * * rows j>0 -> 1/(x_i + y_j) = 1/(g^-i + g^j) * * Since g is primitive and has multiplicative order 255, * g^-i != g^j for the supported indexes where i + j < 255. * * 4. Place a row filled with 1 at the top of the matrix, transforming it * into an Extended Cauchy Matrix. This preserves the nonsingularity * of all square submatrices. * * This first row reproduces the RAID-5 coefficients. * * 5. Adjust each row with a factor so that the first column contains all 1s. * This transformation also preserves the nonsingularity property, * maintaining the MDS code property. * * Concrete example in GF(256) for k=6, m=4: * * First, create a 3x6 Cauchy matrix using x_i = g^-i and y_0 = 0, y_j = g^j for j>0. * For RAID_MODE_CAUCHY_RAID, where g=2, this gives: * * x = { 1, 142, 71, 173, 216, 108 } * y = { 0, 2, 4 } * * The Cauchy matrix is: * * 1 2 4 8 16 32 * 244 83 78 183 118 47 * 167 39 213 59 153 82 * * Divide row 2 by 244 and row 3 by 167. Then extend it with a row of ones * on top, and it remains MDS. This forms the code for m=4, with RAID-6 as * a subset: * * 1 1 1 1 1 1 * 1 2 4 8 16 32 * 1 245 210 196 154 113 * 1 187 166 215 199 7 * * Extending the same computation to k=251 and m=6 gives: * * This results in the matrix A[row, col] defined as: * * 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01... * 01 02 04 08 10 20 40 80 1d 3a 74 e8 cd 87 13 26 4c 98 2d 5a b4 75... * 01 f5 d2 c4 9a 71 f1 7f fc 87 c1 c6 19 2f 40 55 3d ba 53 04 9c 61... * 01 bb a6 d7 c7 07 ce 82 4a 2f a5 9b b6 60 f1 ad e7 f4 06 d2 df 2e... * 01 97 7f 9c 7c 18 bd a2 58 1a da 74 70 a3 e5 47 29 07 f5 80 23 e9... * 01 2b 3f cf 73 2c d6 ed cb 74 15 78 8a c1 17 c9 89 68 21 ab 76 3b... * * This matrix supports six levels of parity, one for each row, for up to 251 * data disks, one for each column. All the 377,342,351,231 square submatrices * are nonsingular, verified also by brute-force testing. * * The RAID/g=2 construction can be extended beyond six parity levels by * adding rows and removing one data column for each additional row, while * preserving the field-size constraint i + j < 255. The current * implementation remains limited to RAID_PARITY_MAX = 6. * (See mktables.c for more details on how the matrix is generated.) * * The downside is that generic multiplications are required to compute the * parity, not only fast multiplications by 2^n or 2^-n, which negatively * affect performance. However, optimized parallel multiplications using SSSE3 * or AVX2 instructions [1][5] make this approach competitive with the triple * parity computation using power coefficients. * * Another advantage of the RAID/g=2 Cauchy matrix is that the first two rows can * replicate the RAID5 and RAID6 approach, resulting in a compatible extension. * SSSE3 or AVX2 instructions are only needed for triple parity or beyond. * * Parity computation is as follows: * * P = sum(Di) * Q = sum(A[1,i] * Di) * R = sum(A[2,i] * Di) * S = sum(A[3,i] * Di) * T = sum(A[4,i] * Di) * U = sum(A[5,i] * Di) for 0 <= i < N * * For RAID_MODE_CAUCHY_RAID, A[1,i]=2^i. For RAID_MODE_CAUCHY_AES, * A[1,i]=q[i] from the G23 sequence. * * Recovery from six disk failures at indices x, y, z, h, v, w (0 <= x < y < z < h < v < w < N) * involves computing parity of the remaining N-6 disks: * * Pa = sum(Di) * Qa = sum(A[1,i] * Di) * Ra = sum(A[2,i] * Di) * Sa = sum(A[3,i] * Di) * Ta = sum(A[4,i] * Di) * Ua = sum(A[5,i] * Di) for i != x,y,z,h,v,w * * Defining: * * Pd = Pa + P * Qd = Qa + Q * Rd = Ra + R * Sd = Sa + S * Td = Ta + T * Ud = Ua + U * * yields: * * Pd = Dx + Dy + Dz + Dh + Dv + Dw * Qd = A[1,x] * Dx + A[1,y] * Dy + A[1,z] * Dz + A[1,h] * Dh + A[1,v] * Dv + A[1,w] * Dw * Rd = A[2,x] * Dx + A[2,y] * Dy + A[2,z] * Dz + A[2,h] * Dh + A[2,v] * Dv + A[2,w] * Dw * Sd = A[3,x] * Dx + A[3,y] * Dy + A[3,z] * Dz + A[3,h] * Dh + A[3,v] * Dv + A[3,w] * Dw * Td = A[4,x] * Dx + A[4,y] * Dy + A[4,z] * Dz + A[4,h] * Dh + A[4,v] * Dv + A[4,w] * Dw * Ud = A[5,x] * Dx + A[5,y] * Dy + A[5,z] * Dz + A[5,h] * Dh + A[5,v] * Dv + A[5,w] * Dw * * This linear system is always solvable since the coefficient matrix is * nonsingular due to the properties of A[]. * * Example performance on a Intel Core i7-10700 CPU @ 2.90GHz, stripe 256 KiB, 8 data disks: * * int8 int32 int64 sse2 sse2e ssse3 ssse3e avx2 avx2e * gen1 24995 46562 62683 77069 * gen2 8243 15879 25788 27292 44907 * genz 5509 10588 14485 13920 25175 * gen3 1190 13201 14307 27268 * gen4 914 10061 11006 21328 * gen5 746 7970 8828 17040 * gen6 630 6627 7516 14584 * * Values are in MiB/s of data processed by a single thread. * Results can be reproduced using "raid/test/speedtest.c". * * Triple parity with power coefficients "1, 2, 2^-1" is slightly faster than * the Cauchy matrix computation if SSSE3 or AVX2 are available: * * int8 int32 int64 sse2 ssse3 avx2 * genz 2337 2874 10920 18944 * * In conclusion, the use of power coefficients, specifically powers * of 1, 2, and 2^-1, is the best option to implement triple parity on CPUs * without SSSE3 and AVX2. * On modern CPUs with SSSE3 or AVX2, the Cauchy matrix is the best * option because it provides a fast and general approach that works for any * number of parities. * * References: * [1] Anvin, "The mathematics of RAID-6", 2004 * [2] MacWilliams, Sloane, "The Theory of Error-Correcting Codes", 1977 * [3] Blomer, "An XOR-Based Erasure-Resilient Coding Scheme", 1995 * [4] Roth, "Introduction to Coding Theory", 2006 * [5] Plank, "Screaming Fast Galois Field Arithmetic Using Intel SIMD Instructions", 2013 */ /** * Buffer filled with 0 used in recovering. */ static void *raid_zero_block; void raid_zero(void *zero) { raid_zero_block = zero; } /** * Inverts the square matrix M of size nxn into V. * * This is not a general matrix inversion because we assume the matrix M * to have all the square submatrix not singular. * We use Gauss elimination to invert. * * @M Matrix to invert with @n rows and @n columns. * @V Destination matrix where the result is put. * @n Number of rows and columns of the matrix. */ void raid_invert(uint8_t *M, uint8_t *V, int n) { int i, j, k; /* set the identity matrix in V */ for (i = 0; i < n; ++i) for (j = 0; j < n; ++j) V[i * n + j] = i == j; /* for each element in the diagonal */ for (k = 0; k < n; ++k) { uint8_t f; /* * The diagonal element cannot be 0 because * we are inverting matrices with all the square * submatrices not singular */ BUG_ON(M[k * n + k] == 0); /* make the diagonal element to be 1 */ f = inv(M[k * n + k]); for (j = 0; j < n; ++j) { M[k * n + j] = mul(f, M[k * n + j]); V[k * n + j] = mul(f, V[k * n + j]); } /* * Make all the elements over and under the diagonal * to be zero */ for (i = 0; i < n; ++i) { if (i == k) continue; f = M[i * n + k]; for (j = 0; j < n; ++j) { M[i * n + j] ^= mul(f, M[k * n + j]); V[i * n + j] ^= mul(f, V[k * n + j]); } } } } /** * Computes the parity without the missing data blocks * and store it in the buffers of such data blocks. * * This is the parity expressed as Pa,Qa,Ra,Sa,Ta,Ua in the equations. */ void raid_delta_gen(int nr, int *id, int *ip, int nd, size_t size, void **v) { void *p[RAID_PARITY_MAX]; void *pa[RAID_PARITY_MAX]; int i, j; int np; void *latest; /* zero buffer must be initialized prior to delta recovery */ BUG_ON(raid_zero_block == 0); /* * Total number of parities we are going to process * they are both the used and the unused ones */ np = ip[nr - 1] + 1; /* latest missing data block */ latest = v[id[nr - 1]]; /* setup pointers for delta computation */ for (i = 0, j = 0; i < np; ++i) { /* keep a copy of the original parity vector */ p[i] = v[nd + i]; if (ip[j] == i) { /* * Set used parities to point to the missing * data blocks. * * The related data blocks are instead set * to point to the "zero" buffer. */ /* * The latest parity to use ends the for loop and * then it cannot happen to process more of them */ BUG_ON(j >= nr); /* buffer for missing data blocks */ pa[j] = v[id[j]]; /* set at zero the missing data blocks */ v[id[j]] = raid_zero_block; /* compute the parity over the missing data blocks */ v[nd + i] = pa[j]; /* check for the next used entry */ ++j; } else { /* * Unused parities are going to be rewritten with * not significant data, because we don't have * functions able to compute only a subset of * parities. * * To avoid this, we reuse parity buffers, * assuming that all the parity functions write * parities in order. * * We assign the unused parity block to the same * block of the latest used parity that we know it * will be written. * * This means that this block will be written * multiple times and only the latest write will * contain the correct data. */ v[nd + i] = latest; } } /* all the parities have to be processed */ BUG_ON(j != nr); /* * Recompute the parity, note that np may be smaller than the * total number of parities available */ raid_gen(nd, np, size, v); /* restore data buffers as before */ for (j = 0; j < nr; ++j) v[id[j]] = pa[j]; /* restore parity buffers as before */ for (i = 0; i < np; ++i) v[nd + i] = p[i]; } /** * Recover failure of one data block for PAR1. * * Starting from the equation: * * Pd = Dx * * and solving we get: * * Dx = Pd */ void raid_rec1of1(int *id, int nd, size_t size, void **v) { void *p; void *pa; /* * For PAR1 we can directly compute the missing block * and we don't need to use the zero buffer */ p = v[nd]; pa = v[id[0]]; /* use the parity as missing data block */ v[id[0]] = p; /* compute the parity over the missing data block */ v[nd] = pa; /* compute */ raid_gen(nd, 1, size, v); /* restore as before */ v[id[0]] = pa; v[nd] = p; } /** * Recover failure of two data blocks for PAR2. * * Starting from the equations: * * Pd = Dx + Dy * Qd = A[1,id[0]] * Dx + A[1,id[1]] * Dy * * and solving we get: * * Qd + A[1,id[0]] * Pd * Dy = ------------------------ * A[1,id[0]] + A[1,id[1]] * * Dx = Dy + Pd * * with conditions: * * A[1,id[0]] + A[1,id[1]] != 0 * * That are always satisfied for any 0<=id[0]