root/drivers/char/ftape/ecc.c

/* */

DEFINITIONS

1. mod255
2. gfadd
3. gfadd_long
4. gfmul
5. gfmul_exp
6. gfmul_exp_long
7. gfdiv
8. gfinv3
9. gfinv2
10. gfmat_mul
11. set_parity
12. compute_syndromes
13. correct_block
14. sanity_check
15. ecc_set_segment_parity
16. ecc_correct_data

   1 /* Yo, Emacs! we're -*- Linux-C -*-
   2  *
   3  *      Copyright (c) 1993 Ning and David Mosberger.
   4  *
   5  * This is based on code originally written by Bas Laarhoven (bas@vimec.nl)
   6  * and David L. Brown, Jr., and incorporates improvements suggested by
   7  * Kai Harrekilde-Petersen.
   8  *
   9  * This program is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU General Public License as
  11  * published by the Free Software Foundation; either version 2, or (at
  12  * your option) any later version.
  13  *
  14  * This program is distributed in the hope that it will be useful, but
  15  * WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU General Public License
  20  * along with this program; see the file COPYING.  If not, write to
  21  * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,
  22  * USA.
  23  *
  24  * $Source: /home/bas/distr/ftape-2.03b/RCS/ecc.c,v $
  25  * $Author: bas $
  26  *
  27  * $Revision: 1.32 $
  28  * $Date: 1995/04/22 07:30:15 $
  29  * $State: Beta $
  30  *
  31  *      This file contains the Reed-Solomon error correction code 
  32  *      for the QIC-40/80 floppy-tape driver for Linux.
  33  */
  34 
  35 #include <linux/ftape.h>
  36 #include <stdio.h>
  37 #include <sys/errno.h>
  38 
  39 #include "tracing.h"
  40 #include "ecc.h"
  41 
  42 /*
  43  * Machines that are big-endian should define macro BIG_ENDIAN.
  44  * Unfortunately, there doesn't appear to be a standard include
  45  * file that works for all OSs.
  46  */
  47 
  48 #if defined(__sparc__) || defined(__hppa)
  49 #define BIG_ENDIAN
  50 #endif                          /* __sparc__ || __hppa */
  51 
  52 #if defined(__mips__)
  53 #error Find a smart way to determine the Endianness of the MIPS CPU
  54 #endif
  55 
  56 #ifdef TEST
  57 
  58 #undef TRACE()
  59 #undef TRACE_()
  60 #undef TRACE()
  61 #undef TRACEi()
  62 #undef TRACElx()
  63 #undef TRACE_FUN()
  64 #undef TRACE_EXIT
  65 #define printk  printf
  66 #define TRACE_FUN( level, name) char __fun[] = name
  67 #define TRACE_EXIT
  68 #define TRACE_(l,m) { if (ftape_ecc_tracing >= (l) && (l) <= TOP_LEVEL) { \
  69     printk( "[%03d] " __FILE__ " (%s) - ", (int)ftape_trace_id++, __fun); \
  70     m;  } }
  71 #define TRACE(l,m) TRACE_(l,printk(m".\n"))
  72 #define TRACEi(l,m,i) TRACE_(l,printk(m" %d.\n",i))
  73 #define TRACElx(l,m,i) TRACE_(l,printk(m" 0x%08lx.\n",i))
  74 
  75 int ftape_ecc_tracing = 1;
  76 unsigned char ftape_trace_id = 0;
  77 
  78 #endif                          /* TEST */
  79 
  80 /*
  81  * Notice: to minimize the potential for confusion, we use r to
  82  *         denote the independent variable of the polynomials
  83  *         in the Galois Field GF(2^8).  We reserve x for polynomials
  84  *         that that have coefficients in GF(2^8).
  85  *         
  86  * The Galois Field in which coefficient arithmetic is performed are
  87  * the polynomials over Z_2 (i.e., 0 and 1) modulo the irreducible
  88  * polynomial f(r), where f(r)=r^8 + r^7 + r^2 + r + 1.  A polynomial
  89  * is represented as a byte with the MSB as the coefficient of r^7 and
  90  * the LSB as the coefficient of r^0.  For example, the binary
  91  * representation of f(x) is 0x187 (of course, this doesn't fit into 8
  92  * bits).  In this field, the polynomial r is a primitive element.
  93  * That is, r^i with i in 0,...,255 enumerates all elements in the
  94  * field.
  95  *
  96  * The generator polynomial for the QIC-80 ECC is
  97  *
  98  *      g(x) = x^3 + r^105*x^2 + r^105*x + 1
  99  *
 100  * which can be factored into:
 101  *
 102  *      g(x) = (x-r^-1)(x-r^0)(x-r^1)
 103  *
 104  * the byte representation of the coefficients are:
 105  *
 106  *      r^105 = 0xc0
 107  *      r^-1  = 0xc3
 108  *      r^0   = 0x01
 109  *      r^1   = 0x02
 110  *
 111  * Notice that r^-1 = r^254 as exponent arithmetic is performed
 112  * modulo 2^8-1 = 255.
 113  *
 114  * For more information on Galois Fields and Reed-Solomon codes,
 115  * refer to any good book.  I found _An Introduction to Error
 116  * Correcting Codes with Applications_ by S. A. Vanstone and
 117  * P. C. van Oorschot to be a good introduction into the former.
 118  * _CODING THEORY: The Essentials_ I found very useful for its
 119  * concise description of Reed-Solomon encoding/decoding.
 120  *
 121  */
 122 
 123 typedef unsigned char Matrix[3][3];
 124 
 125 /*
 126  * gfpow[] is defined such that gfpow[i] returns r^i if
 127  * i is in the range [0..255].
 128  */
 129 static const unsigned char gfpow[] =
 130 {
 131         0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80,
 132         0x87, 0x89, 0x95, 0xad, 0xdd, 0x3d, 0x7a, 0xf4,
 133         0x6f, 0xde, 0x3b, 0x76, 0xec, 0x5f, 0xbe, 0xfb,
 134         0x71, 0xe2, 0x43, 0x86, 0x8b, 0x91, 0xa5, 0xcd,
 135         0x1d, 0x3a, 0x74, 0xe8, 0x57, 0xae, 0xdb, 0x31,
 136         0x62, 0xc4, 0x0f, 0x1e, 0x3c, 0x78, 0xf0, 0x67,
 137         0xce, 0x1b, 0x36, 0x6c, 0xd8, 0x37, 0x6e, 0xdc,
 138         0x3f, 0x7e, 0xfc, 0x7f, 0xfe, 0x7b, 0xf6, 0x6b,
 139         0xd6, 0x2b, 0x56, 0xac, 0xdf, 0x39, 0x72, 0xe4,
 140         0x4f, 0x9e, 0xbb, 0xf1, 0x65, 0xca, 0x13, 0x26,
 141         0x4c, 0x98, 0xb7, 0xe9, 0x55, 0xaa, 0xd3, 0x21,
 142         0x42, 0x84, 0x8f, 0x99, 0xb5, 0xed, 0x5d, 0xba,
 143         0xf3, 0x61, 0xc2, 0x03, 0x06, 0x0c, 0x18, 0x30,
 144         0x60, 0xc0, 0x07, 0x0e, 0x1c, 0x38, 0x70, 0xe0,
 145         0x47, 0x8e, 0x9b, 0xb1, 0xe5, 0x4d, 0x9a, 0xb3,
 146         0xe1, 0x45, 0x8a, 0x93, 0xa1, 0xc5, 0x0d, 0x1a,
 147         0x34, 0x68, 0xd0, 0x27, 0x4e, 0x9c, 0xbf, 0xf9,
 148         0x75, 0xea, 0x53, 0xa6, 0xcb, 0x11, 0x22, 0x44,
 149         0x88, 0x97, 0xa9, 0xd5, 0x2d, 0x5a, 0xb4, 0xef,
 150         0x59, 0xb2, 0xe3, 0x41, 0x82, 0x83, 0x81, 0x85,
 151         0x8d, 0x9d, 0xbd, 0xfd, 0x7d, 0xfa, 0x73, 0xe6,
 152         0x4b, 0x96, 0xab, 0xd1, 0x25, 0x4a, 0x94, 0xaf,
 153         0xd9, 0x35, 0x6a, 0xd4, 0x2f, 0x5e, 0xbc, 0xff,
 154         0x79, 0xf2, 0x63, 0xc6, 0x0b, 0x16, 0x2c, 0x58,
 155         0xb0, 0xe7, 0x49, 0x92, 0xa3, 0xc1, 0x05, 0x0a,
 156         0x14, 0x28, 0x50, 0xa0, 0xc7, 0x09, 0x12, 0x24,
 157         0x48, 0x90, 0xa7, 0xc9, 0x15, 0x2a, 0x54, 0xa8,
 158         0xd7, 0x29, 0x52, 0xa4, 0xcf, 0x19, 0x32, 0x64,
 159         0xc8, 0x17, 0x2e, 0x5c, 0xb8, 0xf7, 0x69, 0xd2,
 160         0x23, 0x46, 0x8c, 0x9f, 0xb9, 0xf5, 0x6d, 0xda,
 161         0x33, 0x66, 0xcc, 0x1f, 0x3e, 0x7c, 0xf8, 0x77,
 162         0xee, 0x5b, 0xb6, 0xeb, 0x51, 0xa2, 0xc3, 0x01
 163 };
 164 
 165 /*
 166  * This is a log table.  That is, gflog[r^i] returns i (modulo f(r)).
 167  * gflog[0] is undefined and the first element is therefore not valid.
 168  */
 169 static const unsigned char gflog[256] =
 170 {
 171         0xff, 0x00, 0x01, 0x63, 0x02, 0xc6, 0x64, 0x6a,
 172         0x03, 0xcd, 0xc7, 0xbc, 0x65, 0x7e, 0x6b, 0x2a,
 173         0x04, 0x8d, 0xce, 0x4e, 0xc8, 0xd4, 0xbd, 0xe1,
 174         0x66, 0xdd, 0x7f, 0x31, 0x6c, 0x20, 0x2b, 0xf3,
 175         0x05, 0x57, 0x8e, 0xe8, 0xcf, 0xac, 0x4f, 0x83,
 176         0xc9, 0xd9, 0xd5, 0x41, 0xbe, 0x94, 0xe2, 0xb4,
 177         0x67, 0x27, 0xde, 0xf0, 0x80, 0xb1, 0x32, 0x35,
 178         0x6d, 0x45, 0x21, 0x12, 0x2c, 0x0d, 0xf4, 0x38,
 179         0x06, 0x9b, 0x58, 0x1a, 0x8f, 0x79, 0xe9, 0x70,
 180         0xd0, 0xc2, 0xad, 0xa8, 0x50, 0x75, 0x84, 0x48,
 181         0xca, 0xfc, 0xda, 0x8a, 0xd6, 0x54, 0x42, 0x24,
 182         0xbf, 0x98, 0x95, 0xf9, 0xe3, 0x5e, 0xb5, 0x15,
 183         0x68, 0x61, 0x28, 0xba, 0xdf, 0x4c, 0xf1, 0x2f,
 184         0x81, 0xe6, 0xb2, 0x3f, 0x33, 0xee, 0x36, 0x10,
 185         0x6e, 0x18, 0x46, 0xa6, 0x22, 0x88, 0x13, 0xf7,
 186         0x2d, 0xb8, 0x0e, 0x3d, 0xf5, 0xa4, 0x39, 0x3b,
 187         0x07, 0x9e, 0x9c, 0x9d, 0x59, 0x9f, 0x1b, 0x08,
 188         0x90, 0x09, 0x7a, 0x1c, 0xea, 0xa0, 0x71, 0x5a,
 189         0xd1, 0x1d, 0xc3, 0x7b, 0xae, 0x0a, 0xa9, 0x91,
 190         0x51, 0x5b, 0x76, 0x72, 0x85, 0xa1, 0x49, 0xeb,
 191         0xcb, 0x7c, 0xfd, 0xc4, 0xdb, 0x1e, 0x8b, 0xd2,
 192         0xd7, 0x92, 0x55, 0xaa, 0x43, 0x0b, 0x25, 0xaf,
 193         0xc0, 0x73, 0x99, 0x77, 0x96, 0x5c, 0xfa, 0x52,
 194         0xe4, 0xec, 0x5f, 0x4a, 0xb6, 0xa2, 0x16, 0x86,
 195         0x69, 0xc5, 0x62, 0xfe, 0x29, 0x7d, 0xbb, 0xcc,
 196         0xe0, 0xd3, 0x4d, 0x8c, 0xf2, 0x1f, 0x30, 0xdc,
 197         0x82, 0xab, 0xe7, 0x56, 0xb3, 0x93, 0x40, 0xd8,
 198         0x34, 0xb0, 0xef, 0x26, 0x37, 0x0c, 0x11, 0x44,
 199         0x6f, 0x78, 0x19, 0x9a, 0x47, 0x74, 0xa7, 0xc1,
 200         0x23, 0x53, 0x89, 0xfb, 0x14, 0x5d, 0xf8, 0x97,
 201         0x2e, 0x4b, 0xb9, 0x60, 0x0f, 0xed, 0x3e, 0xe5,
 202         0xf6, 0x87, 0xa5, 0x17, 0x3a, 0xa3, 0x3c, 0xb7
 203 };
 204 
 205 /*
 206  * This is a multiplication table for the factor
 207  * 0xc0 (i.e., r^105 (modulo f(r)).
 208  * gfmul_c0[f] returns r^105 * f(r) (modulo f(r)).
 209  */
 210 static const unsigned char gfmul_c0[256] =
 211 {
 212         0x00, 0xc0, 0x07, 0xc7, 0x0e, 0xce, 0x09, 0xc9,
 213         0x1c, 0xdc, 0x1b, 0xdb, 0x12, 0xd2, 0x15, 0xd5,
 214         0x38, 0xf8, 0x3f, 0xff, 0x36, 0xf6, 0x31, 0xf1,
 215         0x24, 0xe4, 0x23, 0xe3, 0x2a, 0xea, 0x2d, 0xed,
 216         0x70, 0xb0, 0x77, 0xb7, 0x7e, 0xbe, 0x79, 0xb9,
 217         0x6c, 0xac, 0x6b, 0xab, 0x62, 0xa2, 0x65, 0xa5,
 218         0x48, 0x88, 0x4f, 0x8f, 0x46, 0x86, 0x41, 0x81,
 219         0x54, 0x94, 0x53, 0x93, 0x5a, 0x9a, 0x5d, 0x9d,
 220         0xe0, 0x20, 0xe7, 0x27, 0xee, 0x2e, 0xe9, 0x29,
 221         0xfc, 0x3c, 0xfb, 0x3b, 0xf2, 0x32, 0xf5, 0x35,
 222         0xd8, 0x18, 0xdf, 0x1f, 0xd6, 0x16, 0xd1, 0x11,
 223         0xc4, 0x04, 0xc3, 0x03, 0xca, 0x0a, 0xcd, 0x0d,
 224         0x90, 0x50, 0x97, 0x57, 0x9e, 0x5e, 0x99, 0x59,
 225         0x8c, 0x4c, 0x8b, 0x4b, 0x82, 0x42, 0x85, 0x45,
 226         0xa8, 0x68, 0xaf, 0x6f, 0xa6, 0x66, 0xa1, 0x61,
 227         0xb4, 0x74, 0xb3, 0x73, 0xba, 0x7a, 0xbd, 0x7d,
 228         0x47, 0x87, 0x40, 0x80, 0x49, 0x89, 0x4e, 0x8e,
 229         0x5b, 0x9b, 0x5c, 0x9c, 0x55, 0x95, 0x52, 0x92,
 230         0x7f, 0xbf, 0x78, 0xb8, 0x71, 0xb1, 0x76, 0xb6,
 231         0x63, 0xa3, 0x64, 0xa4, 0x6d, 0xad, 0x6a, 0xaa,
 232         0x37, 0xf7, 0x30, 0xf0, 0x39, 0xf9, 0x3e, 0xfe,
 233         0x2b, 0xeb, 0x2c, 0xec, 0x25, 0xe5, 0x22, 0xe2,
 234         0x0f, 0xcf, 0x08, 0xc8, 0x01, 0xc1, 0x06, 0xc6,
 235         0x13, 0xd3, 0x14, 0xd4, 0x1d, 0xdd, 0x1a, 0xda,
 236         0xa7, 0x67, 0xa0, 0x60, 0xa9, 0x69, 0xae, 0x6e,
 237         0xbb, 0x7b, 0xbc, 0x7c, 0xb5, 0x75, 0xb2, 0x72,
 238         0x9f, 0x5f, 0x98, 0x58, 0x91, 0x51, 0x96, 0x56,
 239         0x83, 0x43, 0x84, 0x44, 0x8d, 0x4d, 0x8a, 0x4a,
 240         0xd7, 0x17, 0xd0, 0x10, 0xd9, 0x19, 0xde, 0x1e,
 241         0xcb, 0x0b, 0xcc, 0x0c, 0xc5, 0x05, 0xc2, 0x02,
 242         0xef, 0x2f, 0xe8, 0x28, 0xe1, 0x21, 0xe6, 0x26,
 243         0xf3, 0x33, 0xf4, 0x34, 0xfd, 0x3d, 0xfa, 0x3a
 244 };
 245 
 246 
 247 /*
 248  * Returns V modulo 255 provided V is in the range -255,-254,...,509.
 249  */
 250 static inline unsigned char mod255(int v)
     /*  */
 251 {
 252         if (v > 0) {
 253                 if (v < 255) {
 254                         return v;
 255                 } else {
 256                         return v - 255;
 257                 }
 258         } else {
 259                 return v + 255;
 260         }
 261 }
 262 
 263 
 264 /*
 265  * Add two numbers in the field.  Addition in this field is
 266  * equivalent to a bit-wise exclusive OR operation---subtraction
 267  * is therefore identical to addition.
 268  */
 269 static inline unsigned char gfadd(unsigned char a, unsigned char b)
     /*  */
 270 {
 271         return a ^ b;
 272 }
 273 
 274 
 275 /*
 276  * Add two vectors of numbers in the field.  Each byte in A and B get
 277  * added individually.
 278  */
 279 static inline unsigned long gfadd_long(unsigned long a, unsigned long b)
     /*  */
 280 {
 281         return a ^ b;
 282 }
 283 
 284 
 285 /*
 286  * Multiply two numbers in the field:
 287  */
 288 static inline unsigned char gfmul(unsigned char a, unsigned char b)
     /*  */
 289 {
 290         if (a && b) {
 291                 return gfpow[mod255(gflog[a] + gflog[b])];
 292         } else {
 293                 return 0;
 294         }
 295 }
 296 
 297 
 298 /*
 299  * Just like gfmul, except we have already looked up the log
 300  * of the second number.
 301  */
 302 static inline unsigned char gfmul_exp(unsigned char a, int b)
     /*  */
 303 {
 304         if (a) {
 305                 return gfpow[mod255(gflog[a] + b)];
 306         } else {
 307                 return 0;
 308         }
 309 }
 310 
 311 
 312 /*
 313  * Just like gfmul_exp, except that A is a vector of numbers.  That is,
 314  * each byte in A gets multiplied by gfpow[mod255(B)].
 315  */
 316 static inline unsigned long gfmul_exp_long(unsigned long a, int b)
     /*  */
 317 {
 318         TRACE_FUN(8, "gfmul_exp_long");
 319         unsigned char t;
 320 
 321         if (sizeof(long) == 4) {
 322                 TRACE_EXIT;
 323                 return
 324                     ((t = a >> 24 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 24) : 0) |
 325                     ((t = a >> 16 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 16) : 0) |
 326                     ((t = a >> 8 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 8) : 0) |
 327                     ((t = a >> 0 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 0) : 0);
 328 #if !defined(linux)
 329         } else if (sizeof(long) == 8) {
 330                 TRACE_EXIT;
 331                 return
 332                     ((t = a >> 56 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 56) : 0) |
 333                     ((t = a >> 48 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 48) : 0) |
 334                     ((t = a >> 40 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 40) : 0) |
 335                     ((t = a >> 32 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 32) : 0) |
 336                     ((t = a >> 24 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 24) : 0) |
 337                     ((t = a >> 16 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 16) : 0) |
 338                     ((t = a >> 8 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 8) : 0) |
 339                     ((t = a >> 0 & 0xff) ? (((unsigned long) gfpow[mod255(gflog[t] + b)]) << 0) : 0);
 340 #endif
 341         } else {
 342                 TRACEx1(1, "Error: size of long is %d bytes", (int) sizeof(long));
 343         }
 344         TRACE_EXIT;
 345         return -1;
 346 }
 347 
 348 
 349 /*
 350  * Divide two numbers in the field.  Returns a/b (modulo f(x)).
 351  */
 352 static inline unsigned char gfdiv(unsigned char a, unsigned char b)
     /*  */
 353 {
 354         TRACE_FUN(8, "gfdiv");
 355         if (!b) {
 356                 TRACE(-1, "Error: division by zero");
 357                 return 0xff;
 358         } else if (a == 0) {
 359                 return 0;
 360         } else {
 361                 return gfpow[mod255(gflog[a] - gflog[b])];
 362         }
 363         TRACE_EXIT;
 364 }
 365 
 366 
 367 /*
 368  * The following functions return the inverse of the matrix of the
 369  * linear system that needs to be solved to determine the error
 370  * magnitudes.  The first deals with matrices of rank 3, while the
 371  * second deals with matrices of rank 2.  The error indices are passed
 372  * in arguments L0,..,L2 (0=first sector, 31=last sector).  The
 373  * error indices must be sorted in ascending order, i.e., L0<L1<L2.
 374  *
 375  * The linear system that needs to be solved for the error
 376  * magnitudes is A * b = s, where s is the known vector of
 377  * syndromes, b is the vector of error magnitudes and A in
 378  * the ORDER=3 case:
 379  *
 380  *    A_3 = {{1/r^L[0], 1/r^L[1], 1/r^L[2]},
 381  *          {        1,        1,        1},
 382  *          {   r^L[0],   r^L[1],   r^L[2]}}
 383  */
 384 static inline int gfinv3(unsigned char l0, unsigned char l1, unsigned char l2, Matrix Ainv)
     /*  */
 385 {
 386         TRACE_FUN(8, "gfinv3");
 387         unsigned char det;
 388         unsigned char t20, t10, t21, t12, t01, t02;
 389         int log_det;
 390 
 391         /* compute some intermediate results: */
 392         t20 = gfpow[l2 - l0];   /* t20 = r^l2/r^l0 */
 393         t10 = gfpow[l1 - l0];   /* t10 = r^l1/r^l0 */
 394         t21 = gfpow[l2 - l1];   /* t21 = r^l2/r^l1 */
 395         t12 = gfpow[l1 - l2 + 255];     /* t12 = r^l1/r^l2 */
 396         t01 = gfpow[l0 - l1 + 255];     /* t01 = r^l0/r^l1 */
 397         t02 = gfpow[l0 - l2 + 255];     /* t02 = r^l0/r^l2 */
 398         /*
 399          * Calculate the determinant of matrix A_3^-1 (sometimes called
 400          * the Vandermonde determinant):
 401          */
 402         det = gfadd(t20, gfadd(t10, gfadd(t21, gfadd(t12, gfadd(t01, t02)))));
 403         if (!det) {
 404                 TRACE(1, "Inversion failed (3 CRC errors, >0 CRC failures)");
 405                 TRACE_EXIT;
 406                 return 0;
 407         }
 408         log_det = 255 - gflog[det];
 409 
 410         /*
 411          * Now, calculate all of the coefficients:
 412          */
 413         Ainv[0][0] = gfmul_exp(gfadd(gfpow[l1], gfpow[l2]), log_det);
 414         Ainv[0][1] = gfmul_exp(gfadd(t21, t12), log_det);
 415         Ainv[0][2] = gfmul_exp(gfadd(gfpow[255 - l1], gfpow[255 - l2]), log_det);
 416 
 417         Ainv[1][0] = gfmul_exp(gfadd(gfpow[l0], gfpow[l2]), log_det);
 418         Ainv[1][1] = gfmul_exp(gfadd(t20, t02), log_det);
 419         Ainv[1][2] = gfmul_exp(gfadd(gfpow[255 - l0], gfpow[255 - l2]), log_det);
 420 
 421         Ainv[2][0] = gfmul_exp(gfadd(gfpow[l0], gfpow[l1]), log_det);
 422         Ainv[2][1] = gfmul_exp(gfadd(t10, t01), log_det);
 423         Ainv[2][2] = gfmul_exp(gfadd(gfpow[255 - l0], gfpow[255 - l1]), log_det);
 424 
 425         TRACE_EXIT;
 426         return 1;
 427 }
 428 
 429 
 430 static inline int gfinv2(unsigned char l0, unsigned char l1, Matrix Ainv)
     /*  */
 431 {
 432         TRACE_FUN(8, "gfinv2");
 433         unsigned char det;
 434         unsigned char t1, t2;
 435         int log_det;
 436 
 437         t1 = gfpow[255 - l0];
 438         t2 = gfpow[255 - l1];
 439         det = gfadd(t1, t2);
 440         if (!det) {
 441                 TRACE(1, "Inversion failed (2 CRC errors, >0 CRC failures)");
 442                 TRACE_EXIT;
 443                 return 0;
 444         }
 445         log_det = 255 - gflog[det];
 446 
 447         /*
 448          * Now, calculate all of the coefficients:
 449          */
 450         Ainv[0][0] = Ainv[1][0] = gfpow[log_det];
 451 
 452         Ainv[0][1] = gfmul_exp(t2, log_det);
 453         Ainv[1][1] = gfmul_exp(t1, log_det);
 454 
 455         TRACE_EXIT;
 456         return 1;
 457 }
 458 
 459 
 460 /*
 461  * Multiply matrix A by vector S and return result in vector B.
 462  * M is assumed to be of order NxN, S and B of order Nx1.
 463  */
 464 static inline void gfmat_mul(int n, Matrix A, unsigned char *s, unsigned char *b)
     /*  */
 465 {
 466         int i, j;
 467         unsigned char dot_prod;
 468 
 469         for (i = 0; i < n; ++i) {
 470                 dot_prod = 0;
 471                 for (j = 0; j < n; ++j) {
 472                         dot_prod = gfadd(dot_prod, gfmul(A[i][j], s[j]));
 473                 }
 474                 b[i] = dot_prod;
 475         }
 476 }
 477 
 478 
 479 
 480 /*
 481  * The Reed Solomon ECC codes are computed over the N-th byte of each
 482  * block, where N=SECTOR_SIZE.  There are up to 29 blocks of data, and
 483  * 3 blocks of ECC.  The blocks are stored contiguously in memory.
 484  * A segment, consequently, is assumed to have at least 4 blocks:
 485  * one or more data blocks plus three ECC blocks.
 486  *
 487  * Notice: In QIC-80 speak, a CRC error is a sector with an incorrect
 488  *         CRC.  A CRC failure is a sector with incorrect data, but
 489  *         a valid CRC.  In the error control literature, the former
 490  *         is usually called "erasure", the latter "error."
 491  */
 492 /*
 493  * Compute the parity bytes for C columns of data, where C is the
 494  * number of bytes that fit into a long integer.  We use a linear
 495  * feed-back register to do this.  The parity bytes P[0], P[STRIDE],
 496  * P[2*STRIDE] are computed such that:
 497  *
 498  *              x^k * p(x) + m(x) = 0 (modulo g(x))
 499  *
 500  * where k = NBLOCKS,
 501  *       p(x) = P[0] + P[STRIDE]*x + P[2*STRIDE]*x^2, and
 502  *       m(x) = sum_{i=0}^k m_i*x^i.
 503  *       m_i  = DATA[i*SECTOR_SIZE]
 504  */
 505 static inline void set_parity(unsigned long *data, int nblocks, unsigned long *p, int stride)
     /*  */
 506 {
 507         TRACE_FUN(8, "set_parity");
 508         unsigned long p0, p1, p2, t1, t2, *end;
 509 
 510         end = data + nblocks * (SECTOR_SIZE / sizeof(long));
 511         p0 = p1 = p2 = 0;
 512         while (data < end) {
 513                 /*
 514                  * The new parity bytes p0_i, p1_i, p2_i are computed from the old
 515                  * values p0_{i-1}, p1_{i-1}, p2_{i-1} recursively as:
 516                  *
 517                  *        p0_i = p1_{i-1} + r^105 * (m_{i-1} - p0_{i-1})
 518                  *        p1_i = p2_{i-1} + r^105 * (m_{i-1} - p0_{i-1})
 519                  *        p2_i =                    (m_{i-1} - p0_{i-1})
 520                  *
 521                  * With the initial condition: p0_0 = p1_0 = p2_0 = 0.
 522                  */
 523                 t1 = gfadd_long(*data, p0);
 524                 /*
 525                  * Multiply each byte in t1 by 0xc0:
 526                  */
 527                 if (sizeof(long) == 4) {
 528                         t2 = ((unsigned long) gfmul_c0[t1 >> 24 & 0xff]) << 24 |
 529                             ((unsigned long) gfmul_c0[t1 >> 16 & 0xff]) << 16 |
 530                             ((unsigned long) gfmul_c0[t1 >> 8 & 0xff]) << 8 |
 531                             ((unsigned long) gfmul_c0[t1 >> 0 & 0xff]) << 0;
 532 #if !defined(linux)
 533                 } else if (sizeof(long) == 8) {
 534                         t2 = ((unsigned long) gfmul_c0[t1 >> 56 & 0xff]) << 56 |
 535                             ((unsigned long) gfmul_c0[t1 >> 48 & 0xff]) << 48 |
 536                             ((unsigned long) gfmul_c0[t1 >> 40 & 0xff]) << 40 |
 537                             ((unsigned long) gfmul_c0[t1 >> 32 & 0xff]) << 32 |
 538                             ((unsigned long) gfmul_c0[t1 >> 24 & 0xff]) << 24 |
 539                             ((unsigned long) gfmul_c0[t1 >> 16 & 0xff]) << 16 |
 540                             ((unsigned long) gfmul_c0[t1 >> 8 & 0xff]) << 8 |
 541                             ((unsigned long) gfmul_c0[t1 >> 0 & 0xff]) << 0;
 542 #endif
 543                 } else {
 544                         TRACEx1(1, "Error: long is of size %d", (int) sizeof(long));
 545                 }
 546                 p0 = gfadd_long(t2, p1);
 547                 p1 = gfadd_long(t2, p2);
 548                 p2 = t1;
 549                 data += SECTOR_SIZE / sizeof(long);
 550         }
 551         *p = p0;
 552         p += stride;
 553         *p = p1;
 554         p += stride;
 555         *p = p2;
 556         TRACE_EXIT;
 557 }
 558 
 559 
 560 /*
 561  * Compute the 3 syndrome values.  DATA should point to the first byte
 562  * of the column for which the syndromes are desired.  The syndromes
 563  * are computed over the first NBLOCKS of rows.  The three bytes will be
 564  * placed in S[0], S[1], and S[2].
 565  *
 566  * S[i] is the value of the "message" polynomial m(x) evaluated at the
 567  * i-th root of the generator polynomial g(x).
 568  *
 569  * As g(x)=(x-r^-1)(x-1)(x-r^1) we evaluate the message polynomial at
 570  * x=r^-1 to get S[0], at x=r^0=1 to get S[1], and at x=r to get S[2].
 571  * This could be done directly and efficiently via the Horner scheme.
 572  * However, it would require multiplication tables for the factors
 573  * r^-1 (0xc3) and r (0x02).  The following scheme does not require
 574  * any multiplication tables beyond what's needed for set_parity()
 575  * anyway and is slightly faster if there are no errors and slightly
 576  * slower if there are errors.  The latter is hopefully the infrequent
 577  * case.
 578  *
 579  * To understand the alternative algorithm, notice that
 580  * set_parity(m, k, p) computes parity bytes such that:
 581  *
 582  *      x^k * p(x) = m(x) (modulo g(x)).
 583  *
 584  * That is, to evaluate m(r^m), where r^m is a root of g(x), we can
 585  * simply evaluate (r^m)^k*p(r^m).  Also, notice that p is 0 if and
 586  * only if s is zero.  That is, if all parity bytes are 0, we know
 587  * there is no error in the data and consequently there is no need to
 588  * compute s(x) at all!  In all other cases, we compute s(x) from p(x)
 589  * by evaluating (r^m)^k*p(r^m) for m=-1, m=0, and m=1.  The p(x)
 590  * polynomial is evaluated via the Horner scheme.
 591  */
 592 static int compute_syndromes(unsigned long *data, int nblocks, unsigned long *s)
     /*  */
 593 {
 594         unsigned long p[3];
 595 
 596         set_parity(data, nblocks, p, 1);
 597         if (p[0] | p[1] | p[2]) {
 598                 /*
 599                  * Some of the checked columns do not have a zero syndrome.  For
 600                  * simplicity, we compute the syndromes for all columns that we
 601                  * have computed the remainders for.
 602                  */
 603                 s[0] = gfmul_exp_long(gfadd_long(p[0], gfmul_exp_long(gfadd_long(p[1],
 604                               gfmul_exp_long(p[2], -1)), -1)), -nblocks);
 605                 s[1] = gfadd_long(gfadd_long(p[2], p[1]), p[0]);
 606                 s[2] = gfmul_exp_long(gfadd_long(p[0], gfmul_exp_long(gfadd_long(p[1],
 607                                  gfmul_exp_long(p[2], 1)), 1)), nblocks);
 608                 return 0;
 609         } else {
 610                 return 1;
 611         }
 612 }
 613 
 614 
 615 /*
 616  * Correct the block in the column pointed to by DATA.  There are NBAD
 617  * CRC errors and their indices are in BAD_LOC[0], up to
 618  * BAD_LOC[NBAD-1].  If NBAD>1, Ainv holds the inverse of the matrix
 619  * of the linear system that needs to be solved to determine the error
 620  * magnitudes.  S[0], S[1], and S[2] are the syndrome values.  If row
 621  * j gets corrected, then bit j will be set in CORRECTION_MAP.
 622  */
 623 static inline int correct_block(unsigned char *data, int nblocks,
     /*  */
 624                                 int nbad, int *bad_loc, Matrix Ainv,
 625                                 unsigned char *s,
 626                                 BAD_SECTOR * correction_map)
 627 {
 628         TRACE_FUN(8, "correct_block");
 629         int ncorrected = 0;
 630         int i;
 631         unsigned char t1, t2;
 632         unsigned char c0, c1, c2;       /* check bytes */
 633         unsigned char error_mag[3], log_error_mag;
 634         unsigned char *dp, l, e;
 635 
 636         switch (nbad) {
 637         case 0:
 638                 /* might have a CRC failure: */
 639                 if (s[0] == 0) {
 640                         /* more than one error */
 641                         TRACE(1, "ECC failed (0 CRC errors, >1 CRC failures)");
 642                         TRACE_EXIT;
 643                         return -1;
 644                 }               /* if */
 645                 t1 = gfdiv(s[1], s[0]);
 646                 if ((bad_loc[nbad++] = gflog[t1]) >= nblocks) {
 647                         TRACE(1, "ECC failed (0 CRC errors, >1 CRC failures): ");
 648                         TRACEi(1, "attempt to correct data at ", bad_loc[0]);
 649                         TRACE_EXIT;
 650                         return -1;
 651                 }
 652                 error_mag[0] = s[1];
 653                 break;
 654         case 1:
 655                 t1 = gfadd(gfmul_exp(s[1], bad_loc[0]), s[2]);
 656                 t2 = gfadd(gfmul_exp(s[0], bad_loc[0]), s[1]);
 657                 if (t1 == 0 && t2 == 0) {
 658                         /* one erasure, no error: */
 659                         Ainv[0][0] = gfpow[bad_loc[0]];
 660                 } else if (t1 == 0 || t2 == 0) {
 661                         /* one erasure and more than one error: */
 662                         TRACE(1, "ECC failed (1 erasure, >1 error)");
 663                         TRACE_EXIT;
 664                         return -1;
 665                 } else {
 666                         /* one erasure, one error: */
 667                         if ((bad_loc[nbad++] = gflog[gfdiv(t1, t2)]) >= nblocks) {
 668                                 TRACE(1, "ECC failed (1 CRC errors, >1 CRC failures): ");
 669                                 TRACEi(1, "attempt to correct data at ", bad_loc[1]);
 670                                 TRACE_EXIT;
 671                                 return -1;
 672                         }       /* if */
 673                         if (!gfinv2(bad_loc[0], bad_loc[1], Ainv)) {
 674                                 /* inversion failed---must have more than one error */
 675                                 TRACE_EXIT;
 676                                 return -1;
 677                         }
 678                 }
 679                 /*
 680                  *  FALL THROUGH TO ERROR MAGNITUDE COMPUTATION:
 681                  */
 682         case 2:
 683         case 3:
 684                 /* compute error magnitudes: */
 685                 gfmat_mul(nbad, Ainv, s, error_mag);
 686                 break;
 687 
 688         default:
 689                 TRACE(1, "Internal Error: number of CRC errors > 3");
 690                 TRACE_EXIT;
 691                 return -1;
 692         }
 693 
 694         /*
 695          * Perform correction by adding ERROR_MAG[i] to the byte at offset
 696          * BAD_LOC[i].  Also add the value of the computed error polynomial
 697          * to the syndrome values.  If the correction was successful, the
 698          * resulting check bytes should be zero (i.e., the corrected data
 699          * is a valid code word).
 700          */
 701         c0 = s[0];
 702         c1 = s[1];
 703         c2 = s[2];
 704         for (i = 0; i < nbad; ++i) {
 705                 e = error_mag[i];
 706                 if (e) {
 707                         /* correct the byte at offset L by magnitude E: */
 708                         l = bad_loc[i];
 709                         dp = &data[l * SECTOR_SIZE];
 710                         *dp = gfadd(*dp, e);
 711                         *correction_map |= 1 << l;
 712                         ++ncorrected;
 713 
 714                         log_error_mag = gflog[e];
 715                         c0 = gfadd(c0, gfpow[mod255(log_error_mag - l)]);
 716                         c1 = gfadd(c1, e);
 717                         c2 = gfadd(c2, gfpow[mod255(log_error_mag + l)]);
 718                 }
 719         }
 720         if (c0 || c1 || c2) {
 721                 TRACE(1, "ECC self-check failed, too many errors");
 722                 TRACE_EXIT;
 723                 return -1;
 724         }
 725         TRACE_EXIT;
 726         return ncorrected;
 727 }
 728 
 729 
 730 #if defined(ECC_SANITY_CHECK) || defined(ECC_PARANOID)
 731 
 732 /*
 733  * Perform a sanity check on the computed parity bytes:
 734  */
 735 static int sanity_check(unsigned long *data, int nblocks)
     /*  */
 736 {
 737         TRACE_FUN(8, "sanity_check");
 738         unsigned long s[3];
 739 
 740         if (!compute_syndromes(data, nblocks, s)) {
 741                 TRACE(-1, "Internal Error: syndrome self-check failed");
 742                 TRACE_EXIT;
 743                 return 0;
 744         }
 745         TRACE_EXIT;
 746         return 1;
 747 }
 748 
 749 #endif                          /* defined(ECC_SANITY_CHECK) || defined(ECC_PARANOID) */
 750 
 751 
 752 
 753 /*
 754  * Compute the parity for an entire segment of data.
 755  */
 756 int ecc_set_segment_parity(struct memory_segment *mseg)
     /*  */
 757 {
 758         int i;
 759         unsigned char *parity_bytes;
 760 
 761         parity_bytes = &mseg->data[(mseg->blocks - 3) * SECTOR_SIZE];
 762         for (i = 0; i < SECTOR_SIZE; i += sizeof(long)) {
 763                 set_parity((unsigned long *) &mseg->data[i], mseg->blocks - 3,
 764                            (unsigned long *) &parity_bytes[i],
 765                            SECTOR_SIZE / sizeof(long));
 766 #ifdef ECC_PARANOID
 767                 if (!sanity_check((unsigned long *) &mseg->data[i], mseg->blocks)) {
 768                         return -1;
 769                 }
 770 #endif                          /* ECC_PARANOID */
 771         }
 772         return 0;
 773 }
 774 
 775 
 776 /*
 777  * Checks and corrects (if possible) the segment MSEG.  Returns one of
 778  * ECC_OK, ECC_CORRECTED, and ECC_FAILED.
 779  */
 780 int ecc_correct_data(struct memory_segment *mseg)
     /*  */
 781 {
 782         TRACE_FUN(5, "ecc_correct_data");
 783         int col, i, result;
 784         int ncorrected = 0;
 785         int nerasures = 0;      /* # of erasures (CRC errors) */
 786         int erasure_loc[3];     /* erasure locations */
 787         unsigned long ss[3];
 788         unsigned char s[3];
 789         Matrix Ainv;
 790 
 791         mseg->corrected = 0;
 792 
 793         /* find first column that has non-zero syndromes: */
 794         for (col = 0; col < SECTOR_SIZE; col += sizeof(long)) {
 795                 if (!compute_syndromes((unsigned long *) &mseg->data[col],
 796                                        mseg->blocks, ss)) {
 797                         /* something is wrong---have to fix things */
 798                         break;
 799                 }
 800         }
 801         if (col >= SECTOR_SIZE) {
 802                 /* all syndromes are ok, therefore nothing to correct */
 803                 TRACE_EXIT;
 804                 return ECC_OK;
 805         }
 806         /* count the number of CRC errors if there were any: */
 807         if (mseg->read_bad) {
 808                 for (i = 0; i < mseg->blocks; i++) {
 809                         if (BAD_CHECK(mseg->read_bad, i)) {
 810                                 if (nerasures >= 3) {
 811                                         /* this is too much for ECC */
 812                                         TRACE(1, "ECC failed (>3 CRC errors)");
 813                                         TRACE_EXIT;
 814                                         return ECC_FAILED;
 815                                 }       /* if */
 816                                 erasure_loc[nerasures++] = i;
 817                         }
 818                 }
 819         }
 820         /*
 821            * If there are at least 2 CRC errors, determine inverse of matrix
 822            * of linear system to be solved:
 823          */
 824         switch (nerasures) {
 825         case 2:
 826                 if (!gfinv2(erasure_loc[0], erasure_loc[1], Ainv)) {
 827                         TRACE_EXIT;
 828                         return ECC_FAILED;
 829                 }
 830                 break;
 831         case 3:
 832                 if (!gfinv3(erasure_loc[0], erasure_loc[1], erasure_loc[2], Ainv)) {
 833                         TRACE_EXIT;
 834                         return ECC_FAILED;
 835                 }
 836                 break;
 837         default:
 838                 /* this is not an error condition... */
 839                 break;
 840         }
 841 
 842         do {
 843                 for (i = 0; i < sizeof(long); ++i) {
 844                         s[0] = ss[0];
 845                         s[1] = ss[1];
 846                         s[2] = ss[2];
 847                         if (s[0] | s[1] | s[2]) {
 848 #ifdef BIG_ENDIAN
 849                                 result = correct_block(&mseg->data[col + sizeof(long) - 1 - i],
 850                                                        mseg->blocks,
 851                                          nerasures, erasure_loc, Ainv, s,
 852                                                        &mseg->corrected);
 853 #else
 854                                 result = correct_block(&mseg->data[col + i], mseg->blocks,
 855                                          nerasures, erasure_loc, Ainv, s,
 856                                                        &mseg->corrected);
 857 #endif
 858                                 if (result < 0) {
 859                                         TRACE_EXIT;
 860                                         return ECC_FAILED;
 861                                 }
 862                                 ncorrected += result;
 863                         }
 864                         ss[0] >>= 8;
 865                         ss[1] >>= 8;
 866                         ss[2] >>= 8;
 867                 }
 868 
 869 #ifdef ECC_SANITY_CHECK
 870                 if (!sanity_check((unsigned long *) &mseg->data[col], mseg->blocks)) {
 871                         TRACE_EXIT;
 872                         return ECC_FAILED;
 873                 }
 874 #endif                          /* ECC_SANITY_CHECK */
 875 
 876                 /* find next column with non-zero syndromes: */
 877                 while ((col += sizeof(long)) < SECTOR_SIZE) {
 878                         if (!compute_syndromes((unsigned long *) &mseg->data[col],
 879                                                mseg->blocks, ss)) {
 880                                 /* something is wrong---have to fix things */
 881                                 break;
 882                         }
 883                 }
 884         } while (col < SECTOR_SIZE);
 885         if (ncorrected && nerasures == 0) {
 886                 TRACE(2, "block contained error not caught by CRC");
 887         }
 888         TRACEi((ncorrected > 0) ? 4 : 8, "number of corrections:", ncorrected);
 889         TRACE_EXIT;
 890         return ncorrected ? ECC_CORRECTED : ECC_OK;
 891 }
 892 
 893 /*** end of ecc.c ***/

/* */