2 * ReedSolomonDecoder.cpp
5 * Created by Christian Brunschen on 05/05/2008.
6 * Copyright 2008 Google UK. All rights reserved.
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
12 * http://www.apache.org/licenses/LICENSE-2.0
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.
24 #include "ReedSolomonDecoder.h"
26 #include "GF256Poly.h"
27 #include "ReedSolomonException.h"
31 namespace reedsolomon {
32 ReedSolomonDecoder::ReedSolomonDecoder(GF256 &fld) : field(fld) {
35 ReedSolomonDecoder::~ReedSolomonDecoder() {
38 void ReedSolomonDecoder::decode(ArrayRef<int> received, int twoS) {
40 cout << "decode(): received = " << &received << ", array = " <<
41 received.array_ << " (" << received.array_->count_ << ")\n";
43 Ref<GF256Poly> poly(new GF256Poly(field, received));
44 cout << "decoding with poly " << *poly << "\n";
45 ArrayRef<int> syndromeCoefficients(new Array<int>(twoS));
47 cout << "syndromeCoefficients array = " <<
48 syndromeCoefficients.array_ << "\n";
51 for (int i = 0; i < twoS; i++) {
52 int eval = poly->evaluateAt(field.exp(i));
53 syndromeCoefficients[syndromeCoefficients->size() - 1 - i] = eval;
60 cout << "before returning: syndromeCoefficients rc = " <<
61 syndromeCoefficients.array_->count_ << "\n";
67 cout << "syndromeCoefficients rc = " <<
68 syndromeCoefficients.array_->count_ << "\n";
70 Ref<GF256Poly> syndrome(new GF256Poly(field, syndromeCoefficients));
72 cout << "syndrome rc = " << syndrome.object_->count_ <<
73 ", syndromeCoefficients rc = " <<
74 syndromeCoefficients.array_->count_ << "\n";
76 Ref<GF256Poly> monomial(field.buildMonomial(twoS, 1));
78 cout << "monopmial rc = " << monomial.object_->count_ << "\n";
80 vector<Ref<GF256Poly> > sigmaOmega
81 (runEuclideanAlgorithm(monomial, syndrome, twoS));
82 ArrayRef<int> errorLocations = findErrorLocations(sigmaOmega[0]);
84 cout << "errorLocations count: " << errorLocations.array_->count_ << "\n";
86 ArrayRef<int> errorMagitudes = findErrorMagnitudes(sigmaOmega[1],
88 for (unsigned i = 0; i < errorLocations->size(); i++) {
89 int position = received->size() - 1 - field.log(errorLocations[i]);
90 received[position] = GF256::addOrSubtract(received[position],
95 vector<Ref<GF256Poly> > ReedSolomonDecoder::
96 runEuclideanAlgorithm(Ref<GF256Poly> a,
99 // Assume a's degree is >= b's
100 if (a->getDegree() < b->getDegree()) {
101 Ref<GF256Poly> tmp = a;
106 Ref<GF256Poly> rLast(a);
108 Ref<GF256Poly> sLast(field.getOne());
109 Ref<GF256Poly> s(field.getZero());
110 Ref<GF256Poly> tLast(field.getZero());
111 Ref<GF256Poly> t(field.getOne());
113 // Run Euclidean algorithm until r's degree is less than R/2
114 while (r->getDegree() >= R / 2) {
115 Ref<GF256Poly> rLastLast(rLast);
116 Ref<GF256Poly> sLastLast(sLast);
117 Ref<GF256Poly> tLastLast(tLast);
122 // Divide rLastLast by rLast, with quotient q and remainder r
123 if (rLast->isZero()) {
124 // Oops, Euclidean algorithm already terminated?
125 throw new ReedSolomonException("r_{i-1} was zero");
128 Ref<GF256Poly> q(field.getZero());
129 int denominatorLeadingTerm = rLast->getCoefficient(rLast->getDegree());
130 int dltInverse = field.inverse(denominatorLeadingTerm);
131 while (r->getDegree() >= rLast->getDegree() && !r->isZero()) {
132 int degreeDiff = r->getDegree() - rLast->getDegree();
133 int scale = field.multiply(r->getCoefficient(r->getDegree()),
135 q = q->addOrSubtract(field.buildMonomial(degreeDiff, scale));
136 r = r->addOrSubtract(rLast->multiplyByMonomial(degreeDiff, scale));
139 s = q->multiply(sLast)->addOrSubtract(sLastLast);
140 t = q->multiply(tLast)->addOrSubtract(tLastLast);
143 int sigmaTildeAtZero = t->getCoefficient(0);
144 if (sigmaTildeAtZero == 0) {
145 throw new ReedSolomonException("sigmaTilde(0) was zero");
148 int inverse = field.inverse(sigmaTildeAtZero);
149 Ref<GF256Poly> sigma(t->multiply(inverse));
150 Ref<GF256Poly> omega(r->multiply(inverse));
152 cout << "t = " << *t << "\n";
153 cout << "r = " << *r << "\n";
154 cout << "sigma = " << *sigma << "\n";
155 cout << "omega = " << *omega << "\n";
157 vector<Ref<GF256Poly> > result(2);
163 ArrayRef<int> ReedSolomonDecoder::
164 findErrorLocations(Ref<GF256Poly> errorLocator) {
165 // This is a direct application of Chien's search
166 int numErrors = errorLocator->getDegree();
167 if (numErrors == 1) { // shortcut
168 ArrayRef<int> result(1);
169 result[0] = errorLocator->getCoefficient(1);
172 ArrayRef<int> result(numErrors);
174 for (int i = 1; i < 256 && e < numErrors; i++) {
175 // cout << "errorLocator(" << i << ") == " << errorLocator->evaluateAt(i) << "\n";
176 if (errorLocator->evaluateAt(i) == 0) {
177 result[e] = field.inverse(i);
181 if (e != numErrors) {
182 throw new ReedSolomonException
183 ("Error locator degree does not match number of roots");
188 ArrayRef<int> ReedSolomonDecoder::
189 findErrorMagnitudes(Ref<GF256Poly> errorEvaluator,
190 ArrayRef<int> errorLocations) {
191 // This is directly applying Forney's Formula
192 int s = errorLocations.size();
193 ArrayRef<int> result(s);
194 for (int i = 0; i < s; i++) {
195 int xiInverse = field.inverse(errorLocations[i]);
197 for (int j = 0; j < s; j++) {
200 field.multiply(denominator,
202 (1, field.multiply(errorLocations[j], xiInverse)));
205 result[i] = field.multiply(errorEvaluator->evaluateAt(xiInverse),
206 field.inverse(denominator));