lib/bch.c

   1 /*
   2  * Generic binary BCH encoding/decoding library
   3  *
   4  * This program is free software; you can redistribute it and/or modify it
   5  * under the terms of the GNU General Public License version 2 as published by
   6  * the Free Software Foundation.
   7  *
   8  * This program is distributed in the hope that it will be useful, but WITHOUT
   9  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
  11  * more details.
  12  *
  13  * You should have received a copy of the GNU General Public License along with
  14  * this program; if not, write to the Free Software Foundation, Inc., 51
  15  * Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  16  *
  17  * Copyright © 2011 Parrot S.A.
  18  *
  19  * Author: Ivan Djelic <ivan.djelic@parrot.com>
  20  *
  21  * Description:
  22  *
  23  * This library provides runtime configurable encoding/decoding of binary
  24  * Bose-Chaudhuri-Hocquenghem (BCH) codes.
  25  *
  26  * Call bch_init to get a pointer to a newly allocated bch_control structure for
  27  * the given m (Galois field order), t (error correction capability) and
  28  * (optional) primitive polynomial parameters.
  29  *
  30  * Call bch_encode to compute and store ecc parity bytes to a given buffer.
  31  * Call bch_decode to detect and locate errors in received data.
  32  *
  33  * On systems supporting hw BCH features, intermediate results may be provided
  34  * to bch_decode in order to skip certain steps. See bch_decode() documentation
  35  * for details.
  36  *
  37  * Option CONFIG_BCH_CONST_PARAMS can be used to force fixed values of
  38  * parameters m and t; thus allowing extra compiler optimizations and providing
  39  * better (up to 2x) encoding performance. Using this option makes sense when
  40  * (m,t) are fixed and known in advance, e.g. when using BCH error correction
  41  * on a particular NAND flash device.
  42  *
  43  * Algorithmic details:
  44  *
  45  * Encoding is performed by processing 32 input bits in parallel, using 4
  46  * remainder lookup tables.
  47  *
  48  * The final stage of decoding involves the following internal steps:
  49  * a. Syndrome computation
  50  * b. Error locator polynomial computation using Berlekamp-Massey algorithm
  51  * c. Error locator root finding (by far the most expensive step)
  52  *
  53  * In this implementation, step c is not performed using the usual Chien search.
  54  * Instead, an alternative approach described in [1] is used. It consists in
  55  * factoring the error locator polynomial using the Berlekamp Trace algorithm
  56  * (BTA) down to a certain degree (4), after which ad hoc low-degree polynomial
  57  * solving techniques [2] are used. The resulting algorithm, called BTZ, yields
  58  * much better performance than Chien search for usual (m,t) values (typically
  59  * m >= 13, t < 32, see [1]).
  60  *
  61  * [1] B. Biswas, V. Herbert. Efficient root finding of polynomials over fields
  62  * of characteristic 2, in: Western European Workshop on Research in Cryptology
  63  * - WEWoRC 2009, Graz, Austria, LNCS, Springer, July 2009, to appear.
  64  * [2] [Zin96] V.A. Zinoviev. On the solution of equations of degree 10 over
  65  * finite fields GF(2^q). In Rapport de recherche INRIA no 2829, 1996.
  66  */
  67
  68 #include <linux/kernel.h>
  69 #include <linux/errno.h>
  70 #include <linux/init.h>
  71 #include <linux/module.h>
  72 #include <linux/slab.h>
  73 #include <linux/bitops.h>
  74 #include <asm/byteorder.h>
  75 #include <linux/bch.h>
  76
  77 #if defined(CONFIG_BCH_CONST_PARAMS)
  78 #define GF_M(_p)               (CONFIG_BCH_CONST_M)
  79 #define GF_T(_p)               (CONFIG_BCH_CONST_T)
  80 #define GF_N(_p)               ((1 << (CONFIG_BCH_CONST_M))-1)
  81 #define BCH_MAX_M              (CONFIG_BCH_CONST_M)
  82 #define BCH_MAX_T              (CONFIG_BCH_CONST_T)
  83 #else
  84 #define GF_M(_p)               ((_p)->m)
  85 #define GF_T(_p)               ((_p)->t)
  86 #define GF_N(_p)               ((_p)->n)
  87 #define BCH_MAX_M              15 /* 2KB */
  88 #define BCH_MAX_T              64 /* 64 bit correction */
  89 #endif
  90
  91 #define BCH_ECC_WORDS(_p)      DIV_ROUND_UP(GF_M(_p)*GF_T(_p), 32)
  92 #define BCH_ECC_BYTES(_p)      DIV_ROUND_UP(GF_M(_p)*GF_T(_p), 8)
  93
  94 #define BCH_ECC_MAX_WORDS      DIV_ROUND_UP(BCH_MAX_M * BCH_MAX_T, 32)
  95
  96 #ifndef dbg
  97 #define dbg(_fmt, args...)     do {} while (0)
  98 #endif
  99
 100 /*
 101  * represent a polynomial over GF(2^m)
 102  */
 103 struct gf_poly {
 104         unsigned int deg;    /* polynomial degree */
 105         unsigned int c[];   /* polynomial terms */
 106 };
 107
 108 /* given its degree, compute a polynomial size in bytes */
 109 #define GF_POLY_SZ(_d) (sizeof(struct gf_poly)+((_d)+1)*sizeof(unsigned int))
 110
 111 /* polynomial of degree 1 */
 112 struct gf_poly_deg1 {
 113         struct gf_poly poly;
 114         unsigned int   c[2];
 115 };
 116
 117 static u8 swap_bits_table[] = {
 118         0x00, 0x80, 0x40, 0xc0, 0x20, 0xa0, 0x60, 0xe0,
 119         0x10, 0x90, 0x50, 0xd0, 0x30, 0xb0, 0x70, 0xf0,
 120         0x08, 0x88, 0x48, 0xc8, 0x28, 0xa8, 0x68, 0xe8,
 121         0x18, 0x98, 0x58, 0xd8, 0x38, 0xb8, 0x78, 0xf8,
 122         0x04, 0x84, 0x44, 0xc4, 0x24, 0xa4, 0x64, 0xe4,
 123         0x14, 0x94, 0x54, 0xd4, 0x34, 0xb4, 0x74, 0xf4,
 124         0x0c, 0x8c, 0x4c, 0xcc, 0x2c, 0xac, 0x6c, 0xec,
 125         0x1c, 0x9c, 0x5c, 0xdc, 0x3c, 0xbc, 0x7c, 0xfc,
 126         0x02, 0x82, 0x42, 0xc2, 0x22, 0xa2, 0x62, 0xe2,
 127         0x12, 0x92, 0x52, 0xd2, 0x32, 0xb2, 0x72, 0xf2,
 128         0x0a, 0x8a, 0x4a, 0xca, 0x2a, 0xaa, 0x6a, 0xea,
 129         0x1a, 0x9a, 0x5a, 0xda, 0x3a, 0xba, 0x7a, 0xfa,
 130         0x06, 0x86, 0x46, 0xc6, 0x26, 0xa6, 0x66, 0xe6,
 131         0x16, 0x96, 0x56, 0xd6, 0x36, 0xb6, 0x76, 0xf6,
 132         0x0e, 0x8e, 0x4e, 0xce, 0x2e, 0xae, 0x6e, 0xee,
 133         0x1e, 0x9e, 0x5e, 0xde, 0x3e, 0xbe, 0x7e, 0xfe,
 134         0x01, 0x81, 0x41, 0xc1, 0x21, 0xa1, 0x61, 0xe1,
 135         0x11, 0x91, 0x51, 0xd1, 0x31, 0xb1, 0x71, 0xf1,
 136         0x09, 0x89, 0x49, 0xc9, 0x29, 0xa9, 0x69, 0xe9,
 137         0x19, 0x99, 0x59, 0xd9, 0x39, 0xb9, 0x79, 0xf9,
 138         0x05, 0x85, 0x45, 0xc5, 0x25, 0xa5, 0x65, 0xe5,
 139         0x15, 0x95, 0x55, 0xd5, 0x35, 0xb5, 0x75, 0xf5,
 140         0x0d, 0x8d, 0x4d, 0xcd, 0x2d, 0xad, 0x6d, 0xed,
 141         0x1d, 0x9d, 0x5d, 0xdd, 0x3d, 0xbd, 0x7d, 0xfd,
 142         0x03, 0x83, 0x43, 0xc3, 0x23, 0xa3, 0x63, 0xe3,
 143         0x13, 0x93, 0x53, 0xd3, 0x33, 0xb3, 0x73, 0xf3,
 144         0x0b, 0x8b, 0x4b, 0xcb, 0x2b, 0xab, 0x6b, 0xeb,
 145         0x1b, 0x9b, 0x5b, 0xdb, 0x3b, 0xbb, 0x7b, 0xfb,
 146         0x07, 0x87, 0x47, 0xc7, 0x27, 0xa7, 0x67, 0xe7,
 147         0x17, 0x97, 0x57, 0xd7, 0x37, 0xb7, 0x77, 0xf7,
 148         0x0f, 0x8f, 0x4f, 0xcf, 0x2f, 0xaf, 0x6f, 0xef,
 149         0x1f, 0x9f, 0x5f, 0xdf, 0x3f, 0xbf, 0x7f, 0xff,
 150 };
 151
 152 static u8 swap_bits(struct bch_control *bch, u8 in)
 153 {
 154         if (!bch->swap_bits)
 155                 return in;
 156
 157         return swap_bits_table[in];
 158 }
 159
 160 /*
 161  * same as bch_encode(), but process input data one byte at a time
 162  */
 163 static void bch_encode_unaligned(struct bch_control *bch,
 164                                  const unsigned char *data, unsigned int len,
 165                                  uint32_t *ecc)
 166 {
 167         int i;
 168         const uint32_t *p;
 169         const int l = BCH_ECC_WORDS(bch)-1;
 170
 171         while (len--) {
 172                 u8 tmp = swap_bits(bch, *data++);
 173
 174                 p = bch->mod8_tab + (l+1)*(((ecc[0] >> 24)^(tmp)) & 0xff);
 175
 176                 for (i = 0; i < l; i++)
 177                         ecc[i] = ((ecc[i] << 8)|(ecc[i+1] >> 24))^(*p++);
 178
 179                 ecc[l] = (ecc[l] << 8)^(*p);
 180         }
 181 }
 182
 183 /*
 184  * convert ecc bytes to aligned, zero-padded 32-bit ecc words
 185  */
 186 static void load_ecc8(struct bch_control *bch, uint32_t *dst,
 187                       const uint8_t *src)
 188 {
 189         uint8_t pad[4] = {0, 0, 0, 0};
 190         unsigned int i, nwords = BCH_ECC_WORDS(bch)-1;
 191
 192         for (i = 0; i < nwords; i++, src += 4)
 193                 dst[i] = ((u32)swap_bits(bch, src[0]) << 24) |
 194                         ((u32)swap_bits(bch, src[1]) << 16) |
 195                         ((u32)swap_bits(bch, src[2]) << 8) |
 196                         swap_bits(bch, src[3]);
 197
 198         memcpy(pad, src, BCH_ECC_BYTES(bch)-4*nwords);
 199         dst[nwords] = ((u32)swap_bits(bch, pad[0]) << 24) |
 200                 ((u32)swap_bits(bch, pad[1]) << 16) |
 201                 ((u32)swap_bits(bch, pad[2]) << 8) |
 202                 swap_bits(bch, pad[3]);
 203 }
 204
 205 /*
 206  * convert 32-bit ecc words to ecc bytes
 207  */
 208 static void store_ecc8(struct bch_control *bch, uint8_t *dst,
 209                        const uint32_t *src)
 210 {
 211         uint8_t pad[4];
 212         unsigned int i, nwords = BCH_ECC_WORDS(bch)-1;
 213
 214         for (i = 0; i < nwords; i++) {
 215                 *dst++ = swap_bits(bch, src[i] >> 24);
 216                 *dst++ = swap_bits(bch, src[i] >> 16);
 217                 *dst++ = swap_bits(bch, src[i] >> 8);
 218                 *dst++ = swap_bits(bch, src[i]);
 219         }
 220         pad[0] = swap_bits(bch, src[nwords] >> 24);
 221         pad[1] = swap_bits(bch, src[nwords] >> 16);
 222         pad[2] = swap_bits(bch, src[nwords] >> 8);
 223         pad[3] = swap_bits(bch, src[nwords]);
 224         memcpy(dst, pad, BCH_ECC_BYTES(bch)-4*nwords);
 225 }
 226
 227 /**
 228  * bch_encode - calculate BCH ecc parity of data
 229  * @bch:   BCH control structure
 230  * @data:  data to encode
 231  * @len:   data length in bytes
 232  * @ecc:   ecc parity data, must be initialized by caller
 233  *
 234  * The @ecc parity array is used both as input and output parameter, in order to
 235  * allow incremental computations. It should be of the size indicated by member
 236  * @ecc_bytes of @bch, and should be initialized to 0 before the first call.
 237  *
 238  * The exact number of computed ecc parity bits is given by member @ecc_bits of
 239  * @bch; it may be less than m*t for large values of t.
 240  */
 241 void bch_encode(struct bch_control *bch, const uint8_t *data,
 242                 unsigned int len, uint8_t *ecc)
 243 {
 244         const unsigned int l = BCH_ECC_WORDS(bch)-1;
 245         unsigned int i, mlen;
 246         unsigned long m;
 247         uint32_t w, r[BCH_ECC_MAX_WORDS];
 248         const size_t r_bytes = BCH_ECC_WORDS(bch) * sizeof(*r);
 249         const uint32_t * const tab0 = bch->mod8_tab;
 250         const uint32_t * const tab1 = tab0 + 256*(l+1);
 251         const uint32_t * const tab2 = tab1 + 256*(l+1);
 252         const uint32_t * const tab3 = tab2 + 256*(l+1);
 253         const uint32_t *pdata, *p0, *p1, *p2, *p3;
 254
 255         if (WARN_ON(r_bytes > sizeof(r)))
 256                 return;
 257
 258         if (ecc) {
 259                 /* load ecc parity bytes into internal 32-bit buffer */
 260                 load_ecc8(bch, bch->ecc_buf, ecc);
 261         } else {
 262                 memset(bch->ecc_buf, 0, r_bytes);
 263         }
 264
 265         /* process first unaligned data bytes */
 266         m = ((unsigned long)data) & 3;
 267         if (m) {
 268                 mlen = (len < (4-m)) ? len : 4-m;
 269                 bch_encode_unaligned(bch, data, mlen, bch->ecc_buf);
 270                 data += mlen;
 271                 len  -= mlen;
 272         }
 273
 274         /* process 32-bit aligned data words */
 275         pdata = (uint32_t *)data;
 276         mlen  = len/4;
 277         data += 4*mlen;
 278         len  -= 4*mlen;
 279         memcpy(r, bch->ecc_buf, r_bytes);
 280
 281         /*
 282          * split each 32-bit word into 4 polynomials of weight 8 as follows:
 283          *
 284          * 31 ...24  23 ...16  15 ... 8  7 ... 0
 285          * xxxxxxxx  yyyyyyyy  zzzzzzzz  tttttttt
 286          *                               tttttttt  mod g = r0 (precomputed)
 287          *                     zzzzzzzz  00000000  mod g = r1 (precomputed)
 288          *           yyyyyyyy  00000000  00000000  mod g = r2 (precomputed)
 289          * xxxxxxxx  00000000  00000000  00000000  mod g = r3 (precomputed)
 290          * xxxxxxxx  yyyyyyyy  zzzzzzzz  tttttttt  mod g = r0^r1^r2^r3
 291          */
 292         while (mlen--) {
 293                 /* input data is read in big-endian format */
 294                 w = cpu_to_be32(*pdata++);
 295                 if (bch->swap_bits)
 296                         w = (u32)swap_bits(bch, w) |
 297                             ((u32)swap_bits(bch, w >> 8) << 8) |
 298                             ((u32)swap_bits(bch, w >> 16) << 16) |
 299                             ((u32)swap_bits(bch, w >> 24) << 24);
 300                 w ^= r[0];
 301                 p0 = tab0 + (l+1)*((w >>  0) & 0xff);
 302                 p1 = tab1 + (l+1)*((w >>  8) & 0xff);
 303                 p2 = tab2 + (l+1)*((w >> 16) & 0xff);
 304                 p3 = tab3 + (l+1)*((w >> 24) & 0xff);
 305
 306                 for (i = 0; i < l; i++)
 307                         r[i] = r[i+1]^p0[i]^p1[i]^p2[i]^p3[i];
 308
 309                 r[l] = p0[l]^p1[l]^p2[l]^p3[l];
 310         }
 311         memcpy(bch->ecc_buf, r, r_bytes);
 312
 313         /* process last unaligned bytes */
 314         if (len)
 315                 bch_encode_unaligned(bch, data, len, bch->ecc_buf);
 316
 317         /* store ecc parity bytes into original parity buffer */
 318         if (ecc)
 319                 store_ecc8(bch, ecc, bch->ecc_buf);
 320 }
 321 EXPORT_SYMBOL_GPL(bch_encode);
 322
 323 static inline int modulo(struct bch_control *bch, unsigned int v)
 324 {
 325         const unsigned int n = GF_N(bch);
 326         while (v >= n) {
 327                 v -= n;
 328                 v = (v & n) + (v >> GF_M(bch));
 329         }
 330         return v;
 331 }
 332
 333 /*
 334  * shorter and faster modulo function, only works when v < 2N.
 335  */
 336 static inline int mod_s(struct bch_control *bch, unsigned int v)
 337 {
 338         const unsigned int n = GF_N(bch);
 339         return (v < n) ? v : v-n;
 340 }
 341
 342 static inline int deg(unsigned int poly)
 343 {
 344         /* polynomial degree is the most-significant bit index */
 345         return fls(poly)-1;
 346 }
 347
 348 static inline int parity(unsigned int x)
 349 {
 350         /*
 351          * public domain code snippet, lifted from
 352          * http://www-graphics.stanford.edu/~seander/bithacks.html
 353          */
 354         x ^= x >> 1;
 355         x ^= x >> 2;
 356         x = (x & 0x11111111U) * 0x11111111U;
 357         return (x >> 28) & 1;
 358 }
 359
 360 /* Galois field basic operations: multiply, divide, inverse, etc. */
 361
 362 static inline unsigned int gf_mul(struct bch_control *bch, unsigned int a,
 363                                   unsigned int b)
 364 {
 365         return (a && b) ? bch->a_pow_tab[mod_s(bch, bch->a_log_tab[a]+
 366                                                bch->a_log_tab[b])] : 0;
 367 }
 368
 369 static inline unsigned int gf_sqr(struct bch_control *bch, unsigned int a)
 370 {
 371         return a ? bch->a_pow_tab[mod_s(bch, 2*bch->a_log_tab[a])] : 0;
 372 }
 373
 374 static inline unsigned int gf_div(struct bch_control *bch, unsigned int a,
 375                                   unsigned int b)
 376 {
 377         return a ? bch->a_pow_tab[mod_s(bch, bch->a_log_tab[a]+
 378                                         GF_N(bch)-bch->a_log_tab[b])] : 0;
 379 }
 380
 381 static inline unsigned int gf_inv(struct bch_control *bch, unsigned int a)
 382 {
 383         return bch->a_pow_tab[GF_N(bch)-bch->a_log_tab[a]];
 384 }
 385
 386 static inline unsigned int a_pow(struct bch_control *bch, int i)
 387 {
 388         return bch->a_pow_tab[modulo(bch, i)];
 389 }
 390
 391 static inline int a_log(struct bch_control *bch, unsigned int x)
 392 {
 393         return bch->a_log_tab[x];
 394 }
 395
 396 static inline int a_ilog(struct bch_control *bch, unsigned int x)
 397 {
 398         return mod_s(bch, GF_N(bch)-bch->a_log_tab[x]);
 399 }
 400
 401 /*
 402  * compute 2t syndromes of ecc polynomial, i.e. ecc(a^j) for j=1..2t
 403  */
 404 static void compute_syndromes(struct bch_control *bch, uint32_t *ecc,
 405                               unsigned int *syn)
 406 {
 407         int i, j, s;
 408         unsigned int m;
 409         uint32_t poly;
 410         const int t = GF_T(bch);
 411
 412         s = bch->ecc_bits;
 413
 414         /* make sure extra bits in last ecc word are cleared */
 415         m = ((unsigned int)s) & 31;
 416         if (m)
 417                 ecc[s/32] &= ~((1u << (32-m))-1);
 418         memset(syn, 0, 2*t*sizeof(*syn));
 419
 420         /* compute v(a^j) for j=1 .. 2t-1 */
 421         do {
 422                 poly = *ecc++;
 423                 s -= 32;
 424                 while (poly) {
 425                         i = deg(poly);
 426                         for (j = 0; j < 2*t; j += 2)
 427                                 syn[j] ^= a_pow(bch, (j+1)*(i+s));
 428
 429                         poly ^= (1 << i);
 430                 }
 431         } while (s > 0);
 432
 433         /* v(a^(2j)) = v(a^j)^2 */
 434         for (j = 0; j < t; j++)
 435                 syn[2*j+1] = gf_sqr(bch, syn[j]);
 436 }
 437
 438 static void gf_poly_copy(struct gf_poly *dst, struct gf_poly *src)
 439 {
 440         memcpy(dst, src, GF_POLY_SZ(src->deg));
 441 }
 442
 443 static int compute_error_locator_polynomial(struct bch_control *bch,
 444                                             const unsigned int *syn)
 445 {
 446         const unsigned int t = GF_T(bch);
 447         const unsigned int n = GF_N(bch);
 448         unsigned int i, j, tmp, l, pd = 1, d = syn[0];
 449         struct gf_poly *elp = bch->elp;
 450         struct gf_poly *pelp = bch->poly_2t[0];
 451         struct gf_poly *elp_copy = bch->poly_2t[1];
 452         int k, pp = -1;
 453
 454         memset(pelp, 0, GF_POLY_SZ(2*t));
 455         memset(elp, 0, GF_POLY_SZ(2*t));
 456
 457         pelp->deg = 0;
 458         pelp->c[0] = 1;
 459         elp->deg = 0;
 460         elp->c[0] = 1;
 461
 462         /* use simplified binary Berlekamp-Massey algorithm */
 463         for (i = 0; (i < t) && (elp->deg <= t); i++) {
 464                 if (d) {
 465                         k = 2*i-pp;
 466                         gf_poly_copy(elp_copy, elp);
 467                         /* e[i+1](X) = e[i](X)+di*dp^-1*X^2(i-p)*e[p](X) */
 468                         tmp = a_log(bch, d)+n-a_log(bch, pd);
 469                         for (j = 0; j <= pelp->deg; j++) {
 470                                 if (pelp->c[j]) {
 471                                         l = a_log(bch, pelp->c[j]);
 472                                         elp->c[j+k] ^= a_pow(bch, tmp+l);
 473                                 }
 474                         }
 475                         /* compute l[i+1] = max(l[i]->c[l[p]+2*(i-p]) */
 476                         tmp = pelp->deg+k;
 477                         if (tmp > elp->deg) {
 478                                 elp->deg = tmp;
 479                                 gf_poly_copy(pelp, elp_copy);
 480                                 pd = d;
 481                                 pp = 2*i;
 482                         }
 483                 }
 484                 /* di+1 = S(2i+3)+elp[i+1].1*S(2i+2)+...+elp[i+1].lS(2i+3-l) */
 485                 if (i < t-1) {
 486                         d = syn[2*i+2];
 487                         for (j = 1; j <= elp->deg; j++)
 488                                 d ^= gf_mul(bch, elp->c[j], syn[2*i+2-j]);
 489                 }
 490         }
 491         dbg("elp=%s\n", gf_poly_str(elp));
 492         return (elp->deg > t) ? -1 : (int)elp->deg;
 493 }
 494
 495 /*
 496  * solve a m x m linear system in GF(2) with an expected number of solutions,
 497  * and return the number of found solutions
 498  */
 499 static int solve_linear_system(struct bch_control *bch, unsigned int *rows,
 500                                unsigned int *sol, int nsol)
 501 {
 502         const int m = GF_M(bch);
 503         unsigned int tmp, mask;
 504         int rem, c, r, p, k, param[BCH_MAX_M];
 505
 506         k = 0;
 507         mask = 1 << m;
 508
 509         /* Gaussian elimination */
 510         for (c = 0; c < m; c++) {
 511                 rem = 0;
 512                 p = c-k;
 513                 /* find suitable row for elimination */
 514                 for (r = p; r < m; r++) {
 515                         if (rows[r] & mask) {
 516                                 if (r != p) {
 517                                         tmp = rows[r];
 518                                         rows[r] = rows[p];
 519                                         rows[p] = tmp;
 520                                 }
 521                                 rem = r+1;
 522                                 break;
 523                         }
 524                 }
 525                 if (rem) {
 526                         /* perform elimination on remaining rows */
 527                         tmp = rows[p];
 528                         for (r = rem; r < m; r++) {
 529                                 if (rows[r] & mask)
 530                                         rows[r] ^= tmp;
 531                         }
 532                 } else {
 533                         /* elimination not needed, store defective row index */
 534                         param[k++] = c;
 535                 }
 536                 mask >>= 1;
 537         }
 538         /* rewrite system, inserting fake parameter rows */
 539         if (k > 0) {
 540                 p = k;
 541                 for (r = m-1; r >= 0; r--) {
 542                         if ((r > m-1-k) && rows[r])
 543                                 /* system has no solution */
 544                                 return 0;
 545
 546                         rows[r] = (p && (r == param[p-1])) ?
 547                                 p--, 1u << (m-r) : rows[r-p];
 548                 }
 549         }
 550
 551         if (nsol != (1 << k))
 552                 /* unexpected number of solutions */
 553                 return 0;
 554
 555         for (p = 0; p < nsol; p++) {
 556                 /* set parameters for p-th solution */
 557                 for (c = 0; c < k; c++)
 558                         rows[param[c]] = (rows[param[c]] & ~1)|((p >> c) & 1);
 559
 560                 /* compute unique solution */
 561                 tmp = 0;
 562                 for (r = m-1; r >= 0; r--) {
 563                         mask = rows[r] & (tmp|1);
 564                         tmp |= parity(mask) << (m-r);
 565                 }
 566                 sol[p] = tmp >> 1;
 567         }
 568         return nsol;
 569 }
 570
 571 /*
 572  * this function builds and solves a linear system for finding roots of a degree
 573  * 4 affine monic polynomial X^4+aX^2+bX+c over GF(2^m).
 574  */
 575 static int find_affine4_roots(struct bch_control *bch, unsigned int a,
 576                               unsigned int b, unsigned int c,
 577                               unsigned int *roots)
 578 {
 579         int i, j, k;
 580         const int m = GF_M(bch);
 581         unsigned int mask = 0xff, t, rows[16] = {0,};
 582
 583         j = a_log(bch, b);
 584         k = a_log(bch, a);
 585         rows[0] = c;
 586
 587         /* build linear system to solve X^4+aX^2+bX+c = 0 */
 588         for (i = 0; i < m; i++) {
 589                 rows[i+1] = bch->a_pow_tab[4*i]^
 590                         (a ? bch->a_pow_tab[mod_s(bch, k)] : 0)^
 591                         (b ? bch->a_pow_tab[mod_s(bch, j)] : 0);
 592                 j++;
 593                 k += 2;
 594         }
 595         /*
 596          * transpose 16x16 matrix before passing it to linear solver
 597          * warning: this code assumes m < 16
 598          */
 599         for (j = 8; j != 0; j >>= 1, mask ^= (mask << j)) {
 600                 for (k = 0; k < 16; k = (k+j+1) & ~j) {
 601                         t = ((rows[k] >> j)^rows[k+j]) & mask;
 602                         rows[k] ^= (t << j);
 603                         rows[k+j] ^= t;
 604                 }
 605         }
 606         return solve_linear_system(bch, rows, roots, 4);
 607 }
 608
 609 /*
 610  * compute root r of a degree 1 polynomial over GF(2^m) (returned as log(1/r))
 611  */
 612 static int find_poly_deg1_roots(struct bch_control *bch, struct gf_poly *poly,
 613                                 unsigned int *roots)
 614 {
 615         int n = 0;
 616
 617         if (poly->c[0])
 618                 /* poly[X] = bX+c with c!=0, root=c/b */
 619                 roots[n++] = mod_s(bch, GF_N(bch)-bch->a_log_tab[poly->c[0]]+
 620                                    bch->a_log_tab[poly->c[1]]);
 621         return n;
 622 }
 623
 624 /*
 625  * compute roots of a degree 2 polynomial over GF(2^m)
 626  */
 627 static int find_poly_deg2_roots(struct bch_control *bch, struct gf_poly *poly,
 628                                 unsigned int *roots)
 629 {
 630         int n = 0, i, l0, l1, l2;
 631         unsigned int u, v, r;
 632
 633         if (poly->c[0] && poly->c[1]) {
 634
 635                 l0 = bch->a_log_tab[poly->c[0]];
 636                 l1 = bch->a_log_tab[poly->c[1]];
 637                 l2 = bch->a_log_tab[poly->c[2]];
 638
 639                 /* using z=a/bX, transform aX^2+bX+c into z^2+z+u (u=ac/b^2) */
 640                 u = a_pow(bch, l0+l2+2*(GF_N(bch)-l1));
 641                 /*
 642                  * let u = sum(li.a^i) i=0..m-1; then compute r = sum(li.xi):
 643                  * r^2+r = sum(li.(xi^2+xi)) = sum(li.(a^i+Tr(a^i).a^k)) =
 644                  * u + sum(li.Tr(a^i).a^k) = u+a^k.Tr(sum(li.a^i)) = u+a^k.Tr(u)
 645                  * i.e. r and r+1 are roots iff Tr(u)=0
 646                  */
 647                 r = 0;
 648                 v = u;
 649                 while (v) {
 650                         i = deg(v);
 651                         r ^= bch->xi_tab[i];
 652                         v ^= (1 << i);
 653                 }
 654                 /* verify root */
 655                 if ((gf_sqr(bch, r)^r) == u) {
 656                         /* reverse z=a/bX transformation and compute log(1/r) */
 657                         roots[n++] = modulo(bch, 2*GF_N(bch)-l1-
 658                                             bch->a_log_tab[r]+l2);
 659                         roots[n++] = modulo(bch, 2*GF_N(bch)-l1-
 660                                             bch->a_log_tab[r^1]+l2);
 661                 }
 662         }
 663         return n;
 664 }
 665
 666 /*
 667  * compute roots of a degree 3 polynomial over GF(2^m)
 668  */
 669 static int find_poly_deg3_roots(struct bch_control *bch, struct gf_poly *poly,
 670                                 unsigned int *roots)
 671 {
 672         int i, n = 0;
 673         unsigned int a, b, c, a2, b2, c2, e3, tmp[4];
 674
 675         if (poly->c[0]) {
 676                 /* transform polynomial into monic X^3 + a2X^2 + b2X + c2 */
 677                 e3 = poly->c[3];
 678                 c2 = gf_div(bch, poly->c[0], e3);
 679                 b2 = gf_div(bch, poly->c[1], e3);
 680                 a2 = gf_div(bch, poly->c[2], e3);
 681
 682                 /* (X+a2)(X^3+a2X^2+b2X+c2) = X^4+aX^2+bX+c (affine) */
 683                 c = gf_mul(bch, a2, c2);           /* c = a2c2      */
 684                 b = gf_mul(bch, a2, b2)^c2;        /* b = a2b2 + c2 */
 685                 a = gf_sqr(bch, a2)^b2;            /* a = a2^2 + b2 */
 686
 687                 /* find the 4 roots of this affine polynomial */
 688                 if (find_affine4_roots(bch, a, b, c, tmp) == 4) {
 689                         /* remove a2 from final list of roots */
 690                         for (i = 0; i < 4; i++) {
 691                                 if (tmp[i] != a2)
 692                                         roots[n++] = a_ilog(bch, tmp[i]);
 693                         }
 694                 }
 695         }
 696         return n;
 697 }
 698
 699 /*
 700  * compute roots of a degree 4 polynomial over GF(2^m)
 701  */
 702 static int find_poly_deg4_roots(struct bch_control *bch, struct gf_poly *poly,
 703                                 unsigned int *roots)
 704 {
 705         int i, l, n = 0;
 706         unsigned int a, b, c, d, e = 0, f, a2, b2, c2, e4;
 707
 708         if (poly->c[0] == 0)
 709                 return 0;
 710
 711         /* transform polynomial into monic X^4 + aX^3 + bX^2 + cX + d */
 712         e4 = poly->c[4];
 713         d = gf_div(bch, poly->c[0], e4);
 714         c = gf_div(bch, poly->c[1], e4);
 715         b = gf_div(bch, poly->c[2], e4);
 716         a = gf_div(bch, poly->c[3], e4);
 717
 718         /* use Y=1/X transformation to get an affine polynomial */
 719         if (a) {
 720                 /* first, eliminate cX by using z=X+e with ae^2+c=0 */
 721                 if (c) {
 722                         /* compute e such that e^2 = c/a */
 723                         f = gf_div(bch, c, a);
 724                         l = a_log(bch, f);
 725                         l += (l & 1) ? GF_N(bch) : 0;
 726                         e = a_pow(bch, l/2);
 727                         /*
 728                          * use transformation z=X+e:
 729                          * z^4+e^4 + a(z^3+ez^2+e^2z+e^3) + b(z^2+e^2) +cz+ce+d
 730                          * z^4 + az^3 + (ae+b)z^2 + (ae^2+c)z+e^4+be^2+ae^3+ce+d
 731                          * z^4 + az^3 + (ae+b)z^2 + e^4+be^2+d
 732                          * z^4 + az^3 +     b'z^2 + d'
 733                          */
 734                         d = a_pow(bch, 2*l)^gf_mul(bch, b, f)^d;
 735                         b = gf_mul(bch, a, e)^b;
 736                 }
 737                 /* now, use Y=1/X to get Y^4 + b/dY^2 + a/dY + 1/d */
 738                 if (d == 0)
 739                         /* assume all roots have multiplicity 1 */
 740                         return 0;
 741
 742                 c2 = gf_inv(bch, d);
 743                 b2 = gf_div(bch, a, d);
 744                 a2 = gf_div(bch, b, d);
 745         } else {
 746                 /* polynomial is already affine */
 747                 c2 = d;
 748                 b2 = c;
 749                 a2 = b;
 750         }
 751         /* find the 4 roots of this affine polynomial */
 752         if (find_affine4_roots(bch, a2, b2, c2, roots) == 4) {
 753                 for (i = 0; i < 4; i++) {
 754                         /* post-process roots (reverse transformations) */
 755                         f = a ? gf_inv(bch, roots[i]) : roots[i];
 756                         roots[i] = a_ilog(bch, f^e);
 757                 }
 758                 n = 4;
 759         }
 760         return n;
 761 }
 762
 763 /*
 764  * build monic, log-based representation of a polynomial
 765  */
 766 static void gf_poly_logrep(struct bch_control *bch,
 767                            const struct gf_poly *a, int *rep)
 768 {
 769         int i, d = a->deg, l = GF_N(bch)-a_log(bch, a->c[a->deg]);
 770
 771         /* represent 0 values with -1; warning, rep[d] is not set to 1 */
 772         for (i = 0; i < d; i++)
 773                 rep[i] = a->c[i] ? mod_s(bch, a_log(bch, a->c[i])+l) : -1;
 774 }
 775
 776 /*
 777  * compute polynomial Euclidean division remainder in GF(2^m)[X]
 778  */
 779 static void gf_poly_mod(struct bch_control *bch, struct gf_poly *a,
 780                         const struct gf_poly *b, int *rep)
 781 {
 782         int la, p, m;
 783         unsigned int i, j, *c = a->c;
 784         const unsigned int d = b->deg;
 785
 786         if (a->deg < d)
 787                 return;
 788
 789         /* reuse or compute log representation of denominator */
 790         if (!rep) {
 791                 rep = bch->cache;
 792                 gf_poly_logrep(bch, b, rep);
 793         }
 794
 795         for (j = a->deg; j >= d; j--) {
 796                 if (c[j]) {
 797                         la = a_log(bch, c[j]);
 798                         p = j-d;
 799                         for (i = 0; i < d; i++, p++) {
 800                                 m = rep[i];
 801                                 if (m >= 0)
 802                                         c[p] ^= bch->a_pow_tab[mod_s(bch,
 803                                                                      m+la)];
 804                         }
 805                 }
 806         }
 807         a->deg = d-1;
 808         while (!c[a->deg] && a->deg)
 809                 a->deg--;
 810 }
 811
 812 /*
 813  * compute polynomial Euclidean division quotient in GF(2^m)[X]
 814  */
 815 static void gf_poly_div(struct bch_control *bch, struct gf_poly *a,
 816                         const struct gf_poly *b, struct gf_poly *q)
 817 {
 818         if (a->deg >= b->deg) {
 819                 q->deg = a->deg-b->deg;
 820                 /* compute a mod b (modifies a) */
 821                 gf_poly_mod(bch, a, b, NULL);
 822                 /* quotient is stored in upper part of polynomial a */
 823                 memcpy(q->c, &a->c[b->deg], (1+q->deg)*sizeof(unsigned int));
 824         } else {
 825                 q->deg = 0;
 826                 q->c[0] = 0;
 827         }
 828 }
 829
 830 /*
 831  * compute polynomial GCD (Greatest Common Divisor) in GF(2^m)[X]
 832  */
 833 static struct gf_poly *gf_poly_gcd(struct bch_control *bch, struct gf_poly *a,
 834                                    struct gf_poly *b)
 835 {
 836         struct gf_poly *tmp;
 837
 838         dbg("gcd(%s,%s)=", gf_poly_str(a), gf_poly_str(b));
 839
 840         if (a->deg < b->deg) {
 841                 tmp = b;
 842                 b = a;
 843                 a = tmp;
 844         }
 845
 846         while (b->deg > 0) {
 847                 gf_poly_mod(bch, a, b, NULL);
 848                 tmp = b;
 849                 b = a;
 850                 a = tmp;
 851         }
 852
 853         dbg("%s\n", gf_poly_str(a));
 854
 855         return a;
 856 }
 857
 858 /*
 859  * Given a polynomial f and an integer k, compute Tr(a^kX) mod f
 860  * This is used in Berlekamp Trace algorithm for splitting polynomials
 861  */
 862 static void compute_trace_bk_mod(struct bch_control *bch, int k,
 863                                  const struct gf_poly *f, struct gf_poly *z,
 864                                  struct gf_poly *out)
 865 {
 866         const int m = GF_M(bch);
 867         int i, j;
 868
 869         /* z contains z^2j mod f */
 870         z->deg = 1;
 871         z->c[0] = 0;
 872         z->c[1] = bch->a_pow_tab[k];
 873
 874         out->deg = 0;
 875         memset(out, 0, GF_POLY_SZ(f->deg));
 876
 877         /* compute f log representation only once */
 878         gf_poly_logrep(bch, f, bch->cache);
 879
 880         for (i = 0; i < m; i++) {
 881                 /* add a^(k*2^i)(z^(2^i) mod f) and compute (z^(2^i) mod f)^2 */
 882                 for (j = z->deg; j >= 0; j--) {
 883                         out->c[j] ^= z->c[j];
 884                         z->c[2*j] = gf_sqr(bch, z->c[j]);
 885                         z->c[2*j+1] = 0;
 886                 }
 887                 if (z->deg > out->deg)
 888                         out->deg = z->deg;
 889
 890                 if (i < m-1) {
 891                         z->deg *= 2;
 892                         /* z^(2(i+1)) mod f = (z^(2^i) mod f)^2 mod f */
 893                         gf_poly_mod(bch, z, f, bch->cache);
 894                 }
 895         }
 896         while (!out->c[out->deg] && out->deg)
 897                 out->deg--;
 898
 899         dbg("Tr(a^%d.X) mod f = %s\n", k, gf_poly_str(out));
 900 }
 901
 902 /*
 903  * factor a polynomial using Berlekamp Trace algorithm (BTA)
 904  */
 905 static void factor_polynomial(struct bch_control *bch, int k, struct gf_poly *f,
 906                               struct gf_poly **g, struct gf_poly **h)
 907 {
 908         struct gf_poly *f2 = bch->poly_2t[0];
 909         struct gf_poly *q  = bch->poly_2t[1];
 910         struct gf_poly *tk = bch->poly_2t[2];
 911         struct gf_poly *z  = bch->poly_2t[3];
 912         struct gf_poly *gcd;
 913
 914         dbg("factoring %s...\n", gf_poly_str(f));
 915
 916         *g = f;
 917         *h = NULL;
 918
 919         /* tk = Tr(a^k.X) mod f */
 920         compute_trace_bk_mod(bch, k, f, z, tk);
 921
 922         if (tk->deg > 0) {
 923                 /* compute g = gcd(f, tk) (destructive operation) */
 924                 gf_poly_copy(f2, f);
 925                 gcd = gf_poly_gcd(bch, f2, tk);
 926                 if (gcd->deg < f->deg) {
 927                         /* compute h=f/gcd(f,tk); this will modify f and q */
 928                         gf_poly_div(bch, f, gcd, q);
 929                         /* store g and h in-place (clobbering f) */
 930                         *h = &((struct gf_poly_deg1 *)f)[gcd->deg].poly;
 931                         gf_poly_copy(*g, gcd);
 932                         gf_poly_copy(*h, q);
 933                 }
 934         }
 935 }
 936
 937 /*
 938  * find roots of a polynomial, using BTZ algorithm; see the beginning of this
 939  * file for details
 940  */
 941 static int find_poly_roots(struct bch_control *bch, unsigned int k,
 942                            struct gf_poly *poly, unsigned int *roots)
 943 {
 944         int cnt;
 945         struct gf_poly *f1, *f2;
 946
 947         switch (poly->deg) {
 948                 /* handle low degree polynomials with ad hoc techniques */
 949         case 1:
 950                 cnt = find_poly_deg1_roots(bch, poly, roots);
 951                 break;
 952         case 2:
 953                 cnt = find_poly_deg2_roots(bch, poly, roots);
 954                 break;
 955         case 3:
 956                 cnt = find_poly_deg3_roots(bch, poly, roots);
 957                 break;
 958         case 4:
 959                 cnt = find_poly_deg4_roots(bch, poly, roots);
 960                 break;
 961         default:
 962                 /* factor polynomial using Berlekamp Trace Algorithm (BTA) */
 963                 cnt = 0;
 964                 if (poly->deg && (k <= GF_M(bch))) {
 965                         factor_polynomial(bch, k, poly, &f1, &f2);
 966                         if (f1)
 967                                 cnt += find_poly_roots(bch, k+1, f1, roots);
 968                         if (f2)
 969                                 cnt += find_poly_roots(bch, k+1, f2, roots+cnt);
 970                 }
 971                 break;
 972         }
 973         return cnt;
 974 }
 975
 976 #if defined(USE_CHIEN_SEARCH)
 977 /*
 978  * exhaustive root search (Chien) implementation - not used, included only for
 979  * reference/comparison tests
 980  */
 981 static int chien_search(struct bch_control *bch, unsigned int len,
 982                         struct gf_poly *p, unsigned int *roots)
 983 {
 984         int m;
 985         unsigned int i, j, syn, syn0, count = 0;
 986         const unsigned int k = 8*len+bch->ecc_bits;
 987
 988         /* use a log-based representation of polynomial */
 989         gf_poly_logrep(bch, p, bch->cache);
 990         bch->cache[p->deg] = 0;
 991         syn0 = gf_div(bch, p->c[0], p->c[p->deg]);
 992
 993         for (i = GF_N(bch)-k+1; i <= GF_N(bch); i++) {
 994                 /* compute elp(a^i) */
 995                 for (j = 1, syn = syn0; j <= p->deg; j++) {
 996                         m = bch->cache[j];
 997                         if (m >= 0)
 998                                 syn ^= a_pow(bch, m+j*i);
 999                 }
1000                 if (syn == 0) {
1001                         roots[count++] = GF_N(bch)-i;
1002                         if (count == p->deg)
1003                                 break;
1004                 }
1005         }
1006         return (count == p->deg) ? count : 0;
1007 }
1008 #define find_poly_roots(_p, _k, _elp, _loc) chien_search(_p, len, _elp, _loc)
1009 #endif /* USE_CHIEN_SEARCH */
1010
1011 /**
1012  * bch_decode - decode received codeword and find bit error locations
1013  * @bch:      BCH control structure
1014  * @data:     received data, ignored if @calc_ecc is provided
1015  * @len:      data length in bytes, must always be provided
1016  * @recv_ecc: received ecc, if NULL then assume it was XORed in @calc_ecc
1017  * @calc_ecc: calculated ecc, if NULL then calc_ecc is computed from @data
1018  * @syn:      hw computed syndrome data (if NULL, syndrome is calculated)
1019  * @errloc:   output array of error locations
1020  *
1021  * Returns:
1022  *  The number of errors found, or -EBADMSG if decoding failed, or -EINVAL if
1023  *  invalid parameters were provided
1024  *
1025  * Depending on the available hw BCH support and the need to compute @calc_ecc
1026  * separately (using bch_encode()), this function should be called with one of
1027  * the following parameter configurations -
1028  *
1029  * by providing @data and @recv_ecc only:
1030  *   bch_decode(@bch, @data, @len, @recv_ecc, NULL, NULL, @errloc)
1031  *
1032  * by providing @recv_ecc and @calc_ecc:
1033  *   bch_decode(@bch, NULL, @len, @recv_ecc, @calc_ecc, NULL, @errloc)
1034  *
1035  * by providing ecc = recv_ecc XOR calc_ecc:
1036  *   bch_decode(@bch, NULL, @len, NULL, ecc, NULL, @errloc)
1037  *
1038  * by providing syndrome results @syn:
1039  *   bch_decode(@bch, NULL, @len, NULL, NULL, @syn, @errloc)
1040  *
1041  * Once bch_decode() has successfully returned with a positive value, error
1042  * locations returned in array @errloc should be interpreted as follows -
1043  *
1044  * if (errloc[n] >= 8*len), then n-th error is located in ecc (no need for
1045  * data correction)
1046  *
1047  * if (errloc[n] < 8*len), then n-th error is located in data and can be
1048  * corrected with statement data[errloc[n]/8] ^= 1 << (errloc[n] % 8);
1049  *
1050  * Note that this function does not perform any data correction by itself, it
1051  * merely indicates error locations.
1052  */
1053 int bch_decode(struct bch_control *bch, const uint8_t *data, unsigned int len,
1054                const uint8_t *recv_ecc, const uint8_t *calc_ecc,
1055                const unsigned int *syn, unsigned int *errloc)
1056 {
1057         const unsigned int ecc_words = BCH_ECC_WORDS(bch);
1058         unsigned int nbits;
1059         int i, err, nroots;
1060         uint32_t sum;
1061
1062         /* sanity check: make sure data length can be handled */
1063         if (8*len > (bch->n-bch->ecc_bits))
1064                 return -EINVAL;
1065
1066         /* if caller does not provide syndromes, compute them */
1067         if (!syn) {
1068                 if (!calc_ecc) {
1069                         /* compute received data ecc into an internal buffer */
1070                         if (!data || !recv_ecc)
1071                                 return -EINVAL;
1072                         bch_encode(bch, data, len, NULL);
1073                 } else {
1074                         /* load provided calculated ecc */
1075                         load_ecc8(bch, bch->ecc_buf, calc_ecc);
1076                 }
1077                 /* load received ecc or assume it was XORed in calc_ecc */
1078                 if (recv_ecc) {
1079                         load_ecc8(bch, bch->ecc_buf2, recv_ecc);
1080                         /* XOR received and calculated ecc */
1081                         for (i = 0, sum = 0; i < (int)ecc_words; i++) {
1082                                 bch->ecc_buf[i] ^= bch->ecc_buf2[i];
1083                                 sum |= bch->ecc_buf[i];
1084                         }
1085                         if (!sum)
1086                                 /* no error found */
1087                                 return 0;
1088                 }
1089                 compute_syndromes(bch, bch->ecc_buf, bch->syn);
1090                 syn = bch->syn;
1091         }
1092
1093         err = compute_error_locator_polynomial(bch, syn);
1094         if (err > 0) {
1095                 nroots = find_poly_roots(bch, 1, bch->elp, errloc);
1096                 if (err != nroots)
1097                         err = -1;
1098         }
1099         if (err > 0) {
1100                 /* post-process raw error locations for easier correction */
1101                 nbits = (len*8)+bch->ecc_bits;
1102                 for (i = 0; i < err; i++) {
1103                         if (errloc[i] >= nbits) {
1104                                 err = -1;
1105                                 break;
1106                         }
1107                         errloc[i] = nbits-1-errloc[i];
1108                         if (!bch->swap_bits)
1109                                 errloc[i] = (errloc[i] & ~7) |
1110                                             (7-(errloc[i] & 7));
1111                 }
1112         }
1113         return (err >= 0) ? err : -EBADMSG;
1114 }
1115 EXPORT_SYMBOL_GPL(bch_decode);
1116
1117 /*
1118  * generate Galois field lookup tables
1119  */
1120 static int build_gf_tables(struct bch_control *bch, unsigned int poly)
1121 {
1122         unsigned int i, x = 1;
1123         const unsigned int k = 1 << deg(poly);
1124
1125         /* primitive polynomial must be of degree m */
1126         if (k != (1u << GF_M(bch)))
1127                 return -1;
1128
1129         for (i = 0; i < GF_N(bch); i++) {
1130                 bch->a_pow_tab[i] = x;
1131                 bch->a_log_tab[x] = i;
1132                 if (i && (x == 1))
1133                         /* polynomial is not primitive (a^i=1 with 0<i<2^m-1) */
1134                         return -1;
1135                 x <<= 1;
1136                 if (x & k)
1137                         x ^= poly;
1138         }
1139         bch->a_pow_tab[GF_N(bch)] = 1;
1140         bch->a_log_tab[0] = 0;
1141
1142         return 0;
1143 }
1144
1145 /*
1146  * compute generator polynomial remainder tables for fast encoding
1147  */
1148 static void build_mod8_tables(struct bch_control *bch, const uint32_t *g)
1149 {
1150         int i, j, b, d;
1151         uint32_t data, hi, lo, *tab;
1152         const int l = BCH_ECC_WORDS(bch);
1153         const int plen = DIV_ROUND_UP(bch->ecc_bits+1, 32);
1154         const int ecclen = DIV_ROUND_UP(bch->ecc_bits, 32);
1155
1156         memset(bch->mod8_tab, 0, 4*256*l*sizeof(*bch->mod8_tab));
1157
1158         for (i = 0; i < 256; i++) {
1159                 /* p(X)=i is a small polynomial of weight <= 8 */
1160                 for (b = 0; b < 4; b++) {
1161                         /* we want to compute (p(X).X^(8*b+deg(g))) mod g(X) */
1162                         tab = bch->mod8_tab + (b*256+i)*l;
1163                         data = i << (8*b);
1164                         while (data) {
1165                                 d = deg(data);
1166                                 /* subtract X^d.g(X) from p(X).X^(8*b+deg(g)) */
1167                                 data ^= g[0] >> (31-d);
1168                                 for (j = 0; j < ecclen; j++) {
1169                                         hi = (d < 31) ? g[j] << (d+1) : 0;
1170                                         lo = (j+1 < plen) ?
1171                                                 g[j+1] >> (31-d) : 0;
1172                                         tab[j] ^= hi|lo;
1173                                 }
1174                         }
1175                 }
1176         }
1177 }
1178
1179 /*
1180  * build a base for factoring degree 2 polynomials
1181  */
1182 static int build_deg2_base(struct bch_control *bch)
1183 {
1184         const int m = GF_M(bch);
1185         int i, j, r;
1186         unsigned int sum, x, y, remaining, ak = 0, xi[BCH_MAX_M];
1187
1188         /* find k s.t. Tr(a^k) = 1 and 0 <= k < m */
1189         for (i = 0; i < m; i++) {
1190                 for (j = 0, sum = 0; j < m; j++)
1191                         sum ^= a_pow(bch, i*(1 << j));
1192
1193                 if (sum) {
1194                         ak = bch->a_pow_tab[i];
1195                         break;
1196                 }
1197         }
1198         /* find xi, i=0..m-1 such that xi^2+xi = a^i+Tr(a^i).a^k */
1199         remaining = m;
1200         memset(xi, 0, sizeof(xi));
1201
1202         for (x = 0; (x <= GF_N(bch)) && remaining; x++) {
1203                 y = gf_sqr(bch, x)^x;
1204                 for (i = 0; i < 2; i++) {
1205                         r = a_log(bch, y);
1206                         if (y && (r < m) && !xi[r]) {
1207                                 bch->xi_tab[r] = x;
1208                                 xi[r] = 1;
1209                                 remaining--;
1210                                 dbg("x%d = %x\n", r, x);
1211                                 break;
1212                         }
1213                         y ^= ak;
1214                 }
1215         }
1216         /* should not happen but check anyway */
1217         return remaining ? -1 : 0;
1218 }
1219
1220 static void *bch_alloc(size_t size, int *err)
1221 {
1222         void *ptr;
1223
1224         ptr = kmalloc(size, GFP_KERNEL);
1225         if (ptr == NULL)
1226                 *err = 1;
1227         return ptr;
1228 }
1229
1230 /*
1231  * compute generator polynomial for given (m,t) parameters.
1232  */
1233 static uint32_t *compute_generator_polynomial(struct bch_control *bch)
1234 {
1235         const unsigned int m = GF_M(bch);
1236         const unsigned int t = GF_T(bch);
1237         int n, err = 0;
1238         unsigned int i, j, nbits, r, word, *roots;
1239         struct gf_poly *g;
1240         uint32_t *genpoly;
1241
1242         g = bch_alloc(GF_POLY_SZ(m*t), &err);
1243         roots = bch_alloc((bch->n+1)*sizeof(*roots), &err);
1244         genpoly = bch_alloc(DIV_ROUND_UP(m*t+1, 32)*sizeof(*genpoly), &err);
1245
1246         if (err) {
1247                 kfree(genpoly);
1248                 genpoly = NULL;
1249                 goto finish;
1250         }
1251
1252         /* enumerate all roots of g(X) */
1253         memset(roots , 0, (bch->n+1)*sizeof(*roots));
1254         for (i = 0; i < t; i++) {
1255                 for (j = 0, r = 2*i+1; j < m; j++) {
1256                         roots[r] = 1;
1257                         r = mod_s(bch, 2*r);
1258                 }
1259         }
1260         /* build generator polynomial g(X) */
1261         g->deg = 0;
1262         g->c[0] = 1;
1263         for (i = 0; i < GF_N(bch); i++) {
1264                 if (roots[i]) {
1265                         /* multiply g(X) by (X+root) */
1266                         r = bch->a_pow_tab[i];
1267                         g->c[g->deg+1] = 1;
1268                         for (j = g->deg; j > 0; j--)
1269                                 g->c[j] = gf_mul(bch, g->c[j], r)^g->c[j-1];
1270
1271                         g->c[0] = gf_mul(bch, g->c[0], r);
1272                         g->deg++;
1273                 }
1274         }
1275         /* store left-justified binary representation of g(X) */
1276         n = g->deg+1;
1277         i = 0;
1278
1279         while (n > 0) {
1280                 nbits = (n > 32) ? 32 : n;
1281                 for (j = 0, word = 0; j < nbits; j++) {
1282                         if (g->c[n-1-j])
1283                                 word |= 1u << (31-j);
1284                 }
1285                 genpoly[i++] = word;
1286                 n -= nbits;
1287         }
1288         bch->ecc_bits = g->deg;
1289
1290 finish:
1291         kfree(g);
1292         kfree(roots);
1293
1294         return genpoly;
1295 }
1296
1297 /**
1298  * bch_init - initialize a BCH encoder/decoder
1299  * @m:          Galois field order, should be in the range 5-15
1300  * @t:          maximum error correction capability, in bits
1301  * @prim_poly:  user-provided primitive polynomial (or 0 to use default)
1302  * @swap_bits:  swap bits within data and syndrome bytes
1303  *
1304  * Returns:
1305  *  a newly allocated BCH control structure if successful, NULL otherwise
1306  *
1307  * This initialization can take some time, as lookup tables are built for fast
1308  * encoding/decoding; make sure not to call this function from a time critical
1309  * path. Usually, bch_init() should be called on module/driver init and
1310  * bch_free() should be called to release memory on exit.
1311  *
1312  * You may provide your own primitive polynomial of degree @m in argument
1313  * @prim_poly, or let bch_init() use its default polynomial.
1314  *
1315  * Once bch_init() has successfully returned a pointer to a newly allocated
1316  * BCH control structure, ecc length in bytes is given by member @ecc_bytes of
1317  * the structure.
1318  */
1319 struct bch_control *bch_init(int m, int t, unsigned int prim_poly,
1320                              bool swap_bits)
1321 {
1322         int err = 0;
1323         unsigned int i, words;
1324         uint32_t *genpoly;
1325         struct bch_control *bch = NULL;
1326
1327         const int min_m = 5;
1328
1329         /* default primitive polynomials */
1330         static const unsigned int prim_poly_tab[] = {
1331                 0x25, 0x43, 0x83, 0x11d, 0x211, 0x409, 0x805, 0x1053, 0x201b,
1332                 0x402b, 0x8003,
1333         };
1334
1335 #if defined(CONFIG_BCH_CONST_PARAMS)
1336         if ((m != (CONFIG_BCH_CONST_M)) || (t != (CONFIG_BCH_CONST_T))) {
1337                 printk(KERN_ERR "bch encoder/decoder was configured to support "
1338                        "parameters m=%d, t=%d only!\n",
1339                        CONFIG_BCH_CONST_M, CONFIG_BCH_CONST_T);
1340                 goto fail;
1341         }
1342 #endif
1343         if ((m < min_m) || (m > BCH_MAX_M))
1344                 /*
1345                  * values of m greater than 15 are not currently supported;
1346                  * supporting m > 15 would require changing table base type
1347                  * (uint16_t) and a small patch in matrix transposition
1348                  */
1349                 goto fail;
1350
1351         if (t > BCH_MAX_T)
1352                 /*
1353                  * we can support larger than 64 bits if necessary, at the
1354                  * cost of higher stack usage.
1355                  */
1356                 goto fail;
1357
1358         /* sanity checks */
1359         if ((t < 1) || (m*t >= ((1 << m)-1)))
1360                 /* invalid t value */
1361                 goto fail;
1362
1363         /* select a primitive polynomial for generating GF(2^m) */
1364         if (prim_poly == 0)
1365                 prim_poly = prim_poly_tab[m-min_m];
1366
1367         bch = kzalloc(sizeof(*bch), GFP_KERNEL);
1368         if (bch == NULL)
1369                 goto fail;
1370
1371         bch->m = m;
1372         bch->t = t;
1373         bch->n = (1 << m)-1;
1374         words  = DIV_ROUND_UP(m*t, 32);
1375         bch->ecc_bytes = DIV_ROUND_UP(m*t, 8);
1376         bch->a_pow_tab = bch_alloc((1+bch->n)*sizeof(*bch->a_pow_tab), &err);
1377         bch->a_log_tab = bch_alloc((1+bch->n)*sizeof(*bch->a_log_tab), &err);
1378         bch->mod8_tab  = bch_alloc(words*1024*sizeof(*bch->mod8_tab), &err);
1379         bch->ecc_buf   = bch_alloc(words*sizeof(*bch->ecc_buf), &err);
1380         bch->ecc_buf2  = bch_alloc(words*sizeof(*bch->ecc_buf2), &err);
1381         bch->xi_tab    = bch_alloc(m*sizeof(*bch->xi_tab), &err);
1382         bch->syn       = bch_alloc(2*t*sizeof(*bch->syn), &err);
1383         bch->cache     = bch_alloc(2*t*sizeof(*bch->cache), &err);
1384         bch->elp       = bch_alloc((t+1)*sizeof(struct gf_poly_deg1), &err);
1385         bch->swap_bits = swap_bits;
1386
1387         for (i = 0; i < ARRAY_SIZE(bch->poly_2t); i++)
1388                 bch->poly_2t[i] = bch_alloc(GF_POLY_SZ(2*t), &err);
1389
1390         if (err)
1391                 goto fail;
1392
1393         err = build_gf_tables(bch, prim_poly);
1394         if (err)
1395                 goto fail;
1396
1397         /* use generator polynomial for computing encoding tables */
1398         genpoly = compute_generator_polynomial(bch);
1399         if (genpoly == NULL)
1400                 goto fail;
1401
1402         build_mod8_tables(bch, genpoly);
1403         kfree(genpoly);
1404
1405         err = build_deg2_base(bch);
1406         if (err)
1407                 goto fail;
1408
1409         return bch;
1410
1411 fail:
1412         bch_free(bch);
1413         return NULL;
1414 }
1415 EXPORT_SYMBOL_GPL(bch_init);
1416
1417 /**
1418  *  bch_free - free the BCH control structure
1419  *  @bch:    BCH control structure to release
1420  */
1421 void bch_free(struct bch_control *bch)
1422 {
1423         unsigned int i;
1424
1425         if (bch) {
1426                 kfree(bch->a_pow_tab);
1427                 kfree(bch->a_log_tab);
1428                 kfree(bch->mod8_tab);
1429                 kfree(bch->ecc_buf);
1430                 kfree(bch->ecc_buf2);
1431                 kfree(bch->xi_tab);
1432                 kfree(bch->syn);
1433                 kfree(bch->cache);
1434                 kfree(bch->elp);
1435
1436                 for (i = 0; i < ARRAY_SIZE(bch->poly_2t); i++)
1437                         kfree(bch->poly_2t[i]);
1438
1439                 kfree(bch);
1440         }
1441 }
1442 EXPORT_SYMBOL_GPL(bch_free);
1443
1444 MODULE_LICENSE("GPL");
1445 MODULE_AUTHOR("Ivan Djelic <ivan.djelic@parrot.com>");
1446 MODULE_DESCRIPTION("Binary BCH encoder/decoder");