cpp/core/src/common/reedsolomon/ReedSolomonDecoder.cpp

   1 /*
   2  *  ReedSolomonDecoder.cpp
   3  *  zxing
   4  *
   5  *  Created by Christian Brunschen on 05/05/2008.
   6  *  Copyright 2008 Google UK. All rights reserved.
   7  *
   8  * Licensed under the Apache License, Version 2.0 (the "License");
   9  * you may not use this file except in compliance with the License.
  10  * You may obtain a copy of the License at
  11  *
  12  *      http://www.apache.org/licenses/LICENSE-2.0
  13  *
  14  * Unless required by applicable law or agreed to in writing, software
  15  * distributed under the License is distributed on an "AS IS" BASIS,
  16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  17  * See the License for the specific language governing permissions and
  18  * limitations under the License.
  19  */
  20
  21 #include <iostream>
  22
  23 #include <memory>
  24 #include "ReedSolomonDecoder.h"
  25 #include "GF256.h"
  26 #include "GF256Poly.h"
  27 #include "ReedSolomonException.h"
  28
  29 using namespace std;
  30
  31 namespace reedsolomon {
  32   ReedSolomonDecoder::ReedSolomonDecoder(GF256 &fld) : field(fld) {
  33   }
  34
  35   ReedSolomonDecoder::~ReedSolomonDecoder() {
  36   }
  37
  38   void ReedSolomonDecoder::decode(ArrayRef<int> received, int twoS) {
  39 #ifdef DEBUG
  40     cout << "decode(): received = " << &received << ", array = " <<
  41       received.array_ << " (" << received.array_->count_ << ")\n";
  42 #endif
  43
  44     Ref<GF256Poly> poly(new GF256Poly(field, received));
  45
  46 #ifdef DEBUG
  47     cout << "decoding with poly " << *poly << "\n";
  48 #endif
  49
  50     ArrayRef<int> syndromeCoefficients(new Array<int>(twoS));
  51
  52 #ifdef DEBUG
  53     cout << "syndromeCoefficients array = " <<
  54       syndromeCoefficients.array_ << "\n";
  55 #endif
  56
  57     bool noError = true;
  58     for (int i = 0; i < twoS; i++) {
  59       int eval = poly->evaluateAt(field.exp(i));
  60       syndromeCoefficients[syndromeCoefficients->size() - 1 - i] = eval;
  61       if (eval != 0) {
  62         noError = false;
  63       }
  64     }
  65     if (noError) {
  66
  67 #ifdef DEBUG
  68       cout << "before returning: syndromeCoefficients rc = " <<
  69         syndromeCoefficients.array_->count_ << "\n";
  70 #endif
  71
  72       return;
  73     }
  74
  75 #ifdef DEBUG
  76     cout << "syndromeCoefficients rc = " <<
  77       syndromeCoefficients.array_->count_ << "\n";
  78 #endif
  79     Ref<GF256Poly> syndrome(new GF256Poly(field, syndromeCoefficients));
  80 #ifdef DEBUG
  81     cout << "syndrome rc = " << syndrome.object_->count_ <<
  82       ", syndromeCoefficients rc = " <<
  83       syndromeCoefficients.array_->count_ << "\n";
  84 #endif
  85     Ref<GF256Poly> monomial(field.buildMonomial(twoS, 1));
  86 #ifdef DEBUG
  87     cout << "monopmial rc = " << monomial.object_->count_ << "\n";
  88 #endif
  89     vector<Ref<GF256Poly> > sigmaOmega
  90       (runEuclideanAlgorithm(monomial, syndrome, twoS));
  91     ArrayRef<int> errorLocations = findErrorLocations(sigmaOmega[0]);
  92 #ifdef DEBUG
  93     cout << "errorLocations count: " << errorLocations.array_->count_ << "\n";
  94 #endif
  95     ArrayRef<int> errorMagitudes = findErrorMagnitudes(sigmaOmega[1],
  96                                                        errorLocations);
  97     for (unsigned i = 0; i < errorLocations->size(); i++) {
  98       int position = received->size() - 1 - field.log(errorLocations[i]);
  99       received[position] = GF256::addOrSubtract(received[position],
 100                                                 errorMagitudes[i]);
 101     }
 102   }
 103
 104   vector<Ref<GF256Poly> > ReedSolomonDecoder::
 105   runEuclideanAlgorithm(Ref<GF256Poly> a,
 106                         Ref<GF256Poly> b,
 107                         int R) {
 108     // Assume a's degree is >= b's
 109     if (a->getDegree() < b->getDegree()) {
 110       Ref<GF256Poly> tmp = a;
 111       a = b;
 112       b = tmp;
 113     }
 114
 115     Ref<GF256Poly> rLast(a);
 116     Ref<GF256Poly> r(b);
 117     Ref<GF256Poly> sLast(field.getOne());
 118     Ref<GF256Poly> s(field.getZero());
 119     Ref<GF256Poly> tLast(field.getZero());
 120     Ref<GF256Poly> t(field.getOne());
 121
 122     // Run Euclidean algorithm until r's degree is less than R/2
 123     while (r->getDegree() >= R / 2) {
 124       Ref<GF256Poly> rLastLast(rLast);
 125       Ref<GF256Poly> sLastLast(sLast);
 126       Ref<GF256Poly> tLastLast(tLast);
 127       rLast = r;
 128       sLast = s;
 129       tLast = t;
 130
 131       // Divide rLastLast by rLast, with quotient q and remainder r
 132       if (rLast->isZero()) {
 133         // Oops, Euclidean algorithm already terminated?
 134         throw new ReedSolomonException("r_{i-1} was zero");
 135       }
 136       r = rLastLast;
 137       Ref<GF256Poly> q(field.getZero());
 138       int denominatorLeadingTerm = rLast->getCoefficient(rLast->getDegree());
 139       int dltInverse = field.inverse(denominatorLeadingTerm);
 140       while (r->getDegree() >= rLast->getDegree() && !r->isZero()) {
 141         int degreeDiff = r->getDegree() - rLast->getDegree();
 142         int scale = field.multiply(r->getCoefficient(r->getDegree()),
 143                                    dltInverse);
 144         q = q->addOrSubtract(field.buildMonomial(degreeDiff, scale));
 145         r = r->addOrSubtract(rLast->multiplyByMonomial(degreeDiff, scale));
 146       }
 147
 148       s = q->multiply(sLast)->addOrSubtract(sLastLast);
 149       t = q->multiply(tLast)->addOrSubtract(tLastLast);
 150     }
 151
 152     int sigmaTildeAtZero = t->getCoefficient(0);
 153     if (sigmaTildeAtZero == 0) {
 154       throw new ReedSolomonException("sigmaTilde(0) was zero");
 155     }
 156
 157     int inverse = field.inverse(sigmaTildeAtZero);
 158     Ref<GF256Poly> sigma(t->multiply(inverse));
 159     Ref<GF256Poly> omega(r->multiply(inverse));
 160
 161 #ifdef DEBUG
 162     cout << "t = " << *t << "\n";
 163     cout << "r = " << *r << "\n";
 164     cout << "sigma = " << *sigma << "\n";
 165     cout << "omega = " << *omega << "\n";
 166 #endif
 167
 168     vector<Ref<GF256Poly> > result(2);
 169     result[0] = sigma;
 170     result[1] = omega;
 171     return result;
 172   }
 173
 174   ArrayRef<int> ReedSolomonDecoder::
 175   findErrorLocations(Ref<GF256Poly> errorLocator) {
 176     // This is a direct application of Chien's search
 177     int numErrors = errorLocator->getDegree();
 178     if (numErrors == 1) { // shortcut
 179       ArrayRef<int> result(1);
 180       result[0] = errorLocator->getCoefficient(1);
 181       return result;
 182     }
 183     ArrayRef<int> result(numErrors);
 184     int e = 0;
 185     for (int i = 1; i < 256 && e < numErrors; i++) {
 186       // cout << "errorLocator(" << i << ") == " << errorLocator->evaluateAt(i) << "\n";
 187       if (errorLocator->evaluateAt(i) == 0) {
 188         result[e] = field.inverse(i);
 189         e++;
 190       }
 191     }
 192     if (e != numErrors) {
 193       throw new ReedSolomonException
 194         ("Error locator degree does not match number of roots");
 195     }
 196     return result;
 197   }
 198
 199   ArrayRef<int> ReedSolomonDecoder::
 200   findErrorMagnitudes(Ref<GF256Poly> errorEvaluator,
 201                       ArrayRef<int> errorLocations) {
 202     // This is directly applying Forney's Formula
 203     int s = errorLocations.size();
 204     ArrayRef<int> result(s);
 205     for (int i = 0; i < s; i++) {
 206       int xiInverse = field.inverse(errorLocations[i]);
 207       int denominator = 1;
 208       for (int j = 0; j < s; j++) {
 209         if (i != j) {
 210           denominator =
 211             field.multiply(denominator,
 212                            GF256::addOrSubtract
 213                             (1, field.multiply(errorLocations[j], xiInverse)));
 214         }
 215       }
 216       result[i] = field.multiply(errorEvaluator->evaluateAt(xiInverse),
 217                                  field.inverse(denominator));
 218     }
 219     return result;
 220   }
 221 }