lib/reed_solomon/decode_rs.c

0001 // SPDX-License-Identifier: GPL-2.0
0002 /*
0003  * Generic Reed Solomon encoder / decoder library
0004  *
0005  * Copyright 2002, Phil Karn, KA9Q
0006  * May be used under the terms of the GNU General Public License (GPL)
0007  *
0008  * Adaption to the kernel by Thomas Gleixner (tglx@linutronix.de)
0009  *
0010  * Generic data width independent code which is included by the wrappers.
0011  */
0012 {
0013     struct rs_codec *rs = rsc->codec;
0014     int deg_lambda, el, deg_omega;
0015     int i, j, r, k, pad;
0016     int nn = rs->nn;
0017     int nroots = rs->nroots;
0018     int fcr = rs->fcr;
0019     int prim = rs->prim;
0020     int iprim = rs->iprim;
0021     uint16_t *alpha_to = rs->alpha_to;
0022     uint16_t *index_of = rs->index_of;
0023     uint16_t u, q, tmp, num1, num2, den, discr_r, syn_error;
0024     int count = 0;
0025     int num_corrected;
0026     uint16_t msk = (uint16_t) rs->nn;
0027
0028     /*
0029      * The decoder buffers are in the rs control struct. They are
0030      * arrays sized [nroots + 1]
0031      */
0032     uint16_t *lambda = rsc->buffers + RS_DECODE_LAMBDA * (nroots + 1);
0033     uint16_t *syn = rsc->buffers + RS_DECODE_SYN * (nroots + 1);
0034     uint16_t *b = rsc->buffers + RS_DECODE_B * (nroots + 1);
0035     uint16_t *t = rsc->buffers + RS_DECODE_T * (nroots + 1);
0036     uint16_t *omega = rsc->buffers + RS_DECODE_OMEGA * (nroots + 1);
0037     uint16_t *root = rsc->buffers + RS_DECODE_ROOT * (nroots + 1);
0038     uint16_t *reg = rsc->buffers + RS_DECODE_REG * (nroots + 1);
0039     uint16_t *loc = rsc->buffers + RS_DECODE_LOC * (nroots + 1);
0040
0041     /* Check length parameter for validity */
0042     pad = nn - nroots - len;
0043     BUG_ON(pad < 0 || pad >= nn - nroots);
0044
0045     /* Does the caller provide the syndrome ? */
0046     if (s != NULL) {
0047         for (i = 0; i < nroots; i++) {
0048             /* The syndrome is in index form,
0049              * so nn represents zero
0050              */
0051             if (s[i] != nn)
0052                 goto decode;
0053         }
0054
0055         /* syndrome is zero, no errors to correct  */
0056         return 0;
0057     }
0058
0059     /* form the syndromes; i.e., evaluate data(x) at roots of
0060      * g(x) */
0061     for (i = 0; i < nroots; i++)
0062         syn[i] = (((uint16_t) data[0]) ^ invmsk) & msk;
0063
0064     for (j = 1; j < len; j++) {
0065         for (i = 0; i < nroots; i++) {
0066             if (syn[i] == 0) {
0067                 syn[i] = (((uint16_t) data[j]) ^
0068                       invmsk) & msk;
0069             } else {
0070                 syn[i] = ((((uint16_t) data[j]) ^
0071                        invmsk) & msk) ^
0072                     alpha_to[rs_modnn(rs, index_of[syn[i]] +
0073                                (fcr + i) * prim)];
0074             }
0075         }
0076     }
0077
0078     for (j = 0; j < nroots; j++) {
0079         for (i = 0; i < nroots; i++) {
0080             if (syn[i] == 0) {
0081                 syn[i] = ((uint16_t) par[j]) & msk;
0082             } else {
0083                 syn[i] = (((uint16_t) par[j]) & msk) ^
0084                     alpha_to[rs_modnn(rs, index_of[syn[i]] +
0085                                (fcr+i)*prim)];
0086             }
0087         }
0088     }
0089     s = syn;
0090
0091     /* Convert syndromes to index form, checking for nonzero condition */
0092     syn_error = 0;
0093     for (i = 0; i < nroots; i++) {
0094         syn_error |= s[i];
0095         s[i] = index_of[s[i]];
0096     }
0097
0098     if (!syn_error) {
0099         /* if syndrome is zero, data[] is a codeword and there are no
0100          * errors to correct. So return data[] unmodified
0101          */
0102         return 0;
0103     }
0104
0105  decode:
0106     memset(&lambda[1], 0, nroots * sizeof(lambda[0]));
0107     lambda[0] = 1;
0108
0109     if (no_eras > 0) {
0110         /* Init lambda to be the erasure locator polynomial */
0111         lambda[1] = alpha_to[rs_modnn(rs,
0112                     prim * (nn - 1 - (eras_pos[0] + pad)))];
0113         for (i = 1; i < no_eras; i++) {
0114             u = rs_modnn(rs, prim * (nn - 1 - (eras_pos[i] + pad)));
0115             for (j = i + 1; j > 0; j--) {
0116                 tmp = index_of[lambda[j - 1]];
0117                 if (tmp != nn) {
0118                     lambda[j] ^=
0119                         alpha_to[rs_modnn(rs, u + tmp)];
0120                 }
0121             }
0122         }
0123     }
0124
0125     for (i = 0; i < nroots + 1; i++)
0126         b[i] = index_of[lambda[i]];
0127
0128     /*
0129      * Begin Berlekamp-Massey algorithm to determine error+erasure
0130      * locator polynomial
0131      */
0132     r = no_eras;
0133     el = no_eras;
0134     while (++r <= nroots) { /* r is the step number */
0135         /* Compute discrepancy at the r-th step in poly-form */
0136         discr_r = 0;
0137         for (i = 0; i < r; i++) {
0138             if ((lambda[i] != 0) && (s[r - i - 1] != nn)) {
0139                 discr_r ^=
0140                     alpha_to[rs_modnn(rs,
0141                               index_of[lambda[i]] +
0142                               s[r - i - 1])];
0143             }
0144         }
0145         discr_r = index_of[discr_r];    /* Index form */
0146         if (discr_r == nn) {
0147             /* 2 lines below: B(x) <-- x*B(x) */
0148             memmove (&b[1], b, nroots * sizeof (b[0]));
0149             b[0] = nn;
0150         } else {
0151             /* 7 lines below: T(x) <-- lambda(x)-discr_r*x*b(x) */
0152             t[0] = lambda[0];
0153             for (i = 0; i < nroots; i++) {
0154                 if (b[i] != nn) {
0155                     t[i + 1] = lambda[i + 1] ^
0156                         alpha_to[rs_modnn(rs, discr_r +
0157                                   b[i])];
0158                 } else
0159                     t[i + 1] = lambda[i + 1];
0160             }
0161             if (2 * el <= r + no_eras - 1) {
0162                 el = r + no_eras - el;
0163                 /*
0164                  * 2 lines below: B(x) <-- inv(discr_r) *
0165                  * lambda(x)
0166                  */
0167                 for (i = 0; i <= nroots; i++) {
0168                     b[i] = (lambda[i] == 0) ? nn :
0169                         rs_modnn(rs, index_of[lambda[i]]
0170                              - discr_r + nn);
0171                 }
0172             } else {
0173                 /* 2 lines below: B(x) <-- x*B(x) */
0174                 memmove(&b[1], b, nroots * sizeof(b[0]));
0175                 b[0] = nn;
0176             }
0177             memcpy(lambda, t, (nroots + 1) * sizeof(t[0]));
0178         }
0179     }
0180
0181     /* Convert lambda to index form and compute deg(lambda(x)) */
0182     deg_lambda = 0;
0183     for (i = 0; i < nroots + 1; i++) {
0184         lambda[i] = index_of[lambda[i]];
0185         if (lambda[i] != nn)
0186             deg_lambda = i;
0187     }
0188
0189     if (deg_lambda == 0) {
0190         /*
0191          * deg(lambda) is zero even though the syndrome is non-zero
0192          * => uncorrectable error detected
0193          */
0194         return -EBADMSG;
0195     }
0196
0197     /* Find roots of error+erasure locator polynomial by Chien search */
0198     memcpy(&reg[1], &lambda[1], nroots * sizeof(reg[0]));
0199     count = 0;      /* Number of roots of lambda(x) */
0200     for (i = 1, k = iprim - 1; i <= nn; i++, k = rs_modnn(rs, k + iprim)) {
0201         q = 1;      /* lambda[0] is always 0 */
0202         for (j = deg_lambda; j > 0; j--) {
0203             if (reg[j] != nn) {
0204                 reg[j] = rs_modnn(rs, reg[j] + j);
0205                 q ^= alpha_to[reg[j]];
0206             }
0207         }
0208         if (q != 0)
0209             continue;   /* Not a root */
0210
0211         if (k < pad) {
0212             /* Impossible error location. Uncorrectable error. */
0213             return -EBADMSG;
0214         }
0215
0216         /* store root (index-form) and error location number */
0217         root[count] = i;
0218         loc[count] = k;
0219         /* If we've already found max possible roots,
0220          * abort the search to save time
0221          */
0222         if (++count == deg_lambda)
0223             break;
0224     }
0225     if (deg_lambda != count) {
0226         /*
0227          * deg(lambda) unequal to number of roots => uncorrectable
0228          * error detected
0229          */
0230         return -EBADMSG;
0231     }
0232     /*
0233      * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
0234      * x**nroots). in index form. Also find deg(omega).
0235      */
0236     deg_omega = deg_lambda - 1;
0237     for (i = 0; i <= deg_omega; i++) {
0238         tmp = 0;
0239         for (j = i; j >= 0; j--) {
0240             if ((s[i - j] != nn) && (lambda[j] != nn))
0241                 tmp ^=
0242                     alpha_to[rs_modnn(rs, s[i - j] + lambda[j])];
0243         }
0244         omega[i] = index_of[tmp];
0245     }
0246
0247     /*
0248      * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
0249      * inv(X(l))**(fcr-1) and den = lambda_pr(inv(X(l))) all in poly-form
0250      * Note: we reuse the buffer for b to store the correction pattern
0251      */
0252     num_corrected = 0;
0253     for (j = count - 1; j >= 0; j--) {
0254         num1 = 0;
0255         for (i = deg_omega; i >= 0; i--) {
0256             if (omega[i] != nn)
0257                 num1 ^= alpha_to[rs_modnn(rs, omega[i] +
0258                             i * root[j])];
0259         }
0260
0261         if (num1 == 0) {
0262             /* Nothing to correct at this position */
0263             b[j] = 0;
0264             continue;
0265         }
0266
0267         num2 = alpha_to[rs_modnn(rs, root[j] * (fcr - 1) + nn)];
0268         den = 0;
0269
0270         /* lambda[i+1] for i even is the formal derivative
0271          * lambda_pr of lambda[i] */
0272         for (i = min(deg_lambda, nroots - 1) & ~1; i >= 0; i -= 2) {
0273             if (lambda[i + 1] != nn) {
0274                 den ^= alpha_to[rs_modnn(rs, lambda[i + 1] +
0275                                i * root[j])];
0276             }
0277         }
0278
0279         b[j] = alpha_to[rs_modnn(rs, index_of[num1] +
0280                            index_of[num2] +
0281                            nn - index_of[den])];
0282         num_corrected++;
0283     }
0284
0285     /*
0286      * We compute the syndrome of the 'error' and check that it matches
0287      * the syndrome of the received word
0288      */
0289     for (i = 0; i < nroots; i++) {
0290         tmp = 0;
0291         for (j = 0; j < count; j++) {
0292             if (b[j] == 0)
0293                 continue;
0294
0295             k = (fcr + i) * prim * (nn-loc[j]-1);
0296             tmp ^= alpha_to[rs_modnn(rs, index_of[b[j]] + k)];
0297         }
0298
0299         if (tmp != alpha_to[s[i]])
0300             return -EBADMSG;
0301     }
0302
0303     /*
0304      * Store the error correction pattern, if a
0305      * correction buffer is available
0306      */
0307     if (corr && eras_pos) {
0308         j = 0;
0309         for (i = 0; i < count; i++) {
0310             if (b[i]) {
0311                 corr[j] = b[i];
0312                 eras_pos[j++] = loc[i] - pad;
0313             }
0314         }
0315     } else if (data && par) {
0316         /* Apply error to data and parity */
0317         for (i = 0; i < count; i++) {
0318             if (loc[i] < (nn - nroots))
0319                 data[loc[i] - pad] ^= b[i];
0320             else
0321                 par[loc[i] - pad - len] ^= b[i];
0322         }
0323     }
0324
0325     return  num_corrected;
0326 }