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

/* */