ReedSolomonCodec.php 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468
  1. <?php
  2. declare(strict_types = 1);
  3. namespace BaconQrCode\Common;
  4. use BaconQrCode\Exception\InvalidArgumentException;
  5. use BaconQrCode\Exception\RuntimeException;
  6. use SplFixedArray;
  7. /**
  8. * Reed-Solomon codec for 8-bit characters.
  9. *
  10. * Based on libfec by Phil Karn, KA9Q.
  11. */
  12. final class ReedSolomonCodec
  13. {
  14. /**
  15. * Symbol size in bits.
  16. *
  17. * @var int
  18. */
  19. private $symbolSize;
  20. /**
  21. * Block size in symbols.
  22. *
  23. * @var int
  24. */
  25. private $blockSize;
  26. /**
  27. * First root of RS code generator polynomial, index form.
  28. *
  29. * @var int
  30. */
  31. private $firstRoot;
  32. /**
  33. * Primitive element to generate polynomial roots, index form.
  34. *
  35. * @var int
  36. */
  37. private $primitive;
  38. /**
  39. * Prim-th root of 1, index form.
  40. *
  41. * @var int
  42. */
  43. private $iPrimitive;
  44. /**
  45. * RS code generator polynomial degree (number of roots).
  46. *
  47. * @var int
  48. */
  49. private $numRoots;
  50. /**
  51. * Padding bytes at front of shortened block.
  52. *
  53. * @var int
  54. */
  55. private $padding;
  56. /**
  57. * Log lookup table.
  58. *
  59. * @var SplFixedArray
  60. */
  61. private $alphaTo;
  62. /**
  63. * Anti-Log lookup table.
  64. *
  65. * @var SplFixedArray
  66. */
  67. private $indexOf;
  68. /**
  69. * Generator polynomial.
  70. *
  71. * @var SplFixedArray
  72. */
  73. private $generatorPoly;
  74. /**
  75. * @throws InvalidArgumentException if symbol size ist not between 0 and 8
  76. * @throws InvalidArgumentException if first root is invalid
  77. * @throws InvalidArgumentException if num roots is invalid
  78. * @throws InvalidArgumentException if padding is invalid
  79. * @throws RuntimeException if field generator polynomial is not primitive
  80. */
  81. public function __construct(
  82. int $symbolSize,
  83. int $gfPoly,
  84. int $firstRoot,
  85. int $primitive,
  86. int $numRoots,
  87. int $padding
  88. ) {
  89. if ($symbolSize < 0 || $symbolSize > 8) {
  90. throw new InvalidArgumentException('Symbol size must be between 0 and 8');
  91. }
  92. if ($firstRoot < 0 || $firstRoot >= (1 << $symbolSize)) {
  93. throw new InvalidArgumentException('First root must be between 0 and ' . (1 << $symbolSize));
  94. }
  95. if ($numRoots < 0 || $numRoots >= (1 << $symbolSize)) {
  96. throw new InvalidArgumentException('Num roots must be between 0 and ' . (1 << $symbolSize));
  97. }
  98. if ($padding < 0 || $padding >= ((1 << $symbolSize) - 1 - $numRoots)) {
  99. throw new InvalidArgumentException(
  100. 'Padding must be between 0 and ' . ((1 << $symbolSize) - 1 - $numRoots)
  101. );
  102. }
  103. $this->symbolSize = $symbolSize;
  104. $this->blockSize = (1 << $symbolSize) - 1;
  105. $this->padding = $padding;
  106. $this->alphaTo = SplFixedArray::fromArray(array_fill(0, $this->blockSize + 1, 0), false);
  107. $this->indexOf = SplFixedArray::fromArray(array_fill(0, $this->blockSize + 1, 0), false);
  108. // Generate galous field lookup table
  109. $this->indexOf[0] = $this->blockSize;
  110. $this->alphaTo[$this->blockSize] = 0;
  111. $sr = 1;
  112. for ($i = 0; $i < $this->blockSize; ++$i) {
  113. $this->indexOf[$sr] = $i;
  114. $this->alphaTo[$i] = $sr;
  115. $sr <<= 1;
  116. if ($sr & (1 << $symbolSize)) {
  117. $sr ^= $gfPoly;
  118. }
  119. $sr &= $this->blockSize;
  120. }
  121. if (1 !== $sr) {
  122. throw new RuntimeException('Field generator polynomial is not primitive');
  123. }
  124. // Form RS code generator polynomial from its roots
  125. $this->generatorPoly = SplFixedArray::fromArray(array_fill(0, $numRoots + 1, 0), false);
  126. $this->firstRoot = $firstRoot;
  127. $this->primitive = $primitive;
  128. $this->numRoots = $numRoots;
  129. // Find prim-th root of 1, used in decoding
  130. for ($iPrimitive = 1; ($iPrimitive % $primitive) !== 0; $iPrimitive += $this->blockSize) {
  131. }
  132. $this->iPrimitive = intdiv($iPrimitive, $primitive);
  133. $this->generatorPoly[0] = 1;
  134. for ($i = 0, $root = $firstRoot * $primitive; $i < $numRoots; ++$i, $root += $primitive) {
  135. $this->generatorPoly[$i + 1] = 1;
  136. for ($j = $i; $j > 0; $j--) {
  137. if ($this->generatorPoly[$j] !== 0) {
  138. $this->generatorPoly[$j] = $this->generatorPoly[$j - 1] ^ $this->alphaTo[
  139. $this->modNn($this->indexOf[$this->generatorPoly[$j]] + $root)
  140. ];
  141. } else {
  142. $this->generatorPoly[$j] = $this->generatorPoly[$j - 1];
  143. }
  144. }
  145. $this->generatorPoly[$j] = $this->alphaTo[$this->modNn($this->indexOf[$this->generatorPoly[0]] + $root)];
  146. }
  147. // Convert generator poly to index form for quicker encoding
  148. for ($i = 0; $i <= $numRoots; ++$i) {
  149. $this->generatorPoly[$i] = $this->indexOf[$this->generatorPoly[$i]];
  150. }
  151. }
  152. /**
  153. * Encodes data and writes result back into parity array.
  154. */
  155. public function encode(SplFixedArray $data, SplFixedArray $parity) : void
  156. {
  157. for ($i = 0; $i < $this->numRoots; ++$i) {
  158. $parity[$i] = 0;
  159. }
  160. $iterations = $this->blockSize - $this->numRoots - $this->padding;
  161. for ($i = 0; $i < $iterations; ++$i) {
  162. $feedback = $this->indexOf[$data[$i] ^ $parity[0]];
  163. if ($feedback !== $this->blockSize) {
  164. // Feedback term is non-zero
  165. $feedback = $this->modNn($this->blockSize - $this->generatorPoly[$this->numRoots] + $feedback);
  166. for ($j = 1; $j < $this->numRoots; ++$j) {
  167. $parity[$j] = $parity[$j] ^ $this->alphaTo[
  168. $this->modNn($feedback + $this->generatorPoly[$this->numRoots - $j])
  169. ];
  170. }
  171. }
  172. for ($j = 0; $j < $this->numRoots - 1; ++$j) {
  173. $parity[$j] = $parity[$j + 1];
  174. }
  175. if ($feedback !== $this->blockSize) {
  176. $parity[$this->numRoots - 1] = $this->alphaTo[$this->modNn($feedback + $this->generatorPoly[0])];
  177. } else {
  178. $parity[$this->numRoots - 1] = 0;
  179. }
  180. }
  181. }
  182. /**
  183. * Decodes received data.
  184. */
  185. public function decode(SplFixedArray $data, SplFixedArray $erasures = null) : ?int
  186. {
  187. // This speeds up the initialization a bit.
  188. $numRootsPlusOne = SplFixedArray::fromArray(array_fill(0, $this->numRoots + 1, 0), false);
  189. $numRoots = SplFixedArray::fromArray(array_fill(0, $this->numRoots, 0), false);
  190. $lambda = clone $numRootsPlusOne;
  191. $b = clone $numRootsPlusOne;
  192. $t = clone $numRootsPlusOne;
  193. $omega = clone $numRootsPlusOne;
  194. $root = clone $numRoots;
  195. $loc = clone $numRoots;
  196. $numErasures = (null !== $erasures ? count($erasures) : 0);
  197. // Form the Syndromes; i.e., evaluate data(x) at roots of g(x)
  198. $syndromes = SplFixedArray::fromArray(array_fill(0, $this->numRoots, $data[0]), false);
  199. for ($i = 1; $i < $this->blockSize - $this->padding; ++$i) {
  200. for ($j = 0; $j < $this->numRoots; ++$j) {
  201. if ($syndromes[$j] === 0) {
  202. $syndromes[$j] = $data[$i];
  203. } else {
  204. $syndromes[$j] = $data[$i] ^ $this->alphaTo[
  205. $this->modNn($this->indexOf[$syndromes[$j]] + ($this->firstRoot + $j) * $this->primitive)
  206. ];
  207. }
  208. }
  209. }
  210. // Convert syndromes to index form, checking for nonzero conditions
  211. $syndromeError = 0;
  212. for ($i = 0; $i < $this->numRoots; ++$i) {
  213. $syndromeError |= $syndromes[$i];
  214. $syndromes[$i] = $this->indexOf[$syndromes[$i]];
  215. }
  216. if (! $syndromeError) {
  217. // If syndrome is zero, data[] is a codeword and there are no errors to correct, so return data[]
  218. // unmodified.
  219. return 0;
  220. }
  221. $lambda[0] = 1;
  222. if ($numErasures > 0) {
  223. // Init lambda to be the erasure locator polynomial
  224. $lambda[1] = $this->alphaTo[$this->modNn($this->primitive * ($this->blockSize - 1 - $erasures[0]))];
  225. for ($i = 1; $i < $numErasures; ++$i) {
  226. $u = $this->modNn($this->primitive * ($this->blockSize - 1 - $erasures[$i]));
  227. for ($j = $i + 1; $j > 0; --$j) {
  228. $tmp = $this->indexOf[$lambda[$j - 1]];
  229. if ($tmp !== $this->blockSize) {
  230. $lambda[$j] = $lambda[$j] ^ $this->alphaTo[$this->modNn($u + $tmp)];
  231. }
  232. }
  233. }
  234. }
  235. for ($i = 0; $i <= $this->numRoots; ++$i) {
  236. $b[$i] = $this->indexOf[$lambda[$i]];
  237. }
  238. // Begin Berlekamp-Massey algorithm to determine error+erasure locator polynomial
  239. $r = $numErasures;
  240. $el = $numErasures;
  241. while (++$r <= $this->numRoots) {
  242. // Compute discrepancy at the r-th step in poly form
  243. $discrepancyR = 0;
  244. for ($i = 0; $i < $r; ++$i) {
  245. if ($lambda[$i] !== 0 && $syndromes[$r - $i - 1] !== $this->blockSize) {
  246. $discrepancyR ^= $this->alphaTo[
  247. $this->modNn($this->indexOf[$lambda[$i]] + $syndromes[$r - $i - 1])
  248. ];
  249. }
  250. }
  251. $discrepancyR = $this->indexOf[$discrepancyR];
  252. if ($discrepancyR === $this->blockSize) {
  253. $tmp = $b->toArray();
  254. array_unshift($tmp, $this->blockSize);
  255. array_pop($tmp);
  256. $b = SplFixedArray::fromArray($tmp, false);
  257. continue;
  258. }
  259. $t[0] = $lambda[0];
  260. for ($i = 0; $i < $this->numRoots; ++$i) {
  261. if ($b[$i] !== $this->blockSize) {
  262. $t[$i + 1] = $lambda[$i + 1] ^ $this->alphaTo[$this->modNn($discrepancyR + $b[$i])];
  263. } else {
  264. $t[$i + 1] = $lambda[$i + 1];
  265. }
  266. }
  267. if (2 * $el <= $r + $numErasures - 1) {
  268. $el = $r + $numErasures - $el;
  269. for ($i = 0; $i <= $this->numRoots; ++$i) {
  270. $b[$i] = (
  271. $lambda[$i] === 0
  272. ? $this->blockSize
  273. : $this->modNn($this->indexOf[$lambda[$i]] - $discrepancyR + $this->blockSize)
  274. );
  275. }
  276. } else {
  277. $tmp = $b->toArray();
  278. array_unshift($tmp, $this->blockSize);
  279. array_pop($tmp);
  280. $b = SplFixedArray::fromArray($tmp, false);
  281. }
  282. $lambda = clone $t;
  283. }
  284. // Convert lambda to index form and compute deg(lambda(x))
  285. $degLambda = 0;
  286. for ($i = 0; $i <= $this->numRoots; ++$i) {
  287. $lambda[$i] = $this->indexOf[$lambda[$i]];
  288. if ($lambda[$i] !== $this->blockSize) {
  289. $degLambda = $i;
  290. }
  291. }
  292. // Find roots of the error+erasure locator polynomial by Chien search.
  293. $reg = clone $lambda;
  294. $reg[0] = 0;
  295. $count = 0;
  296. $i = 1;
  297. for ($k = $this->iPrimitive - 1; $i <= $this->blockSize; ++$i, $k = $this->modNn($k + $this->iPrimitive)) {
  298. $q = 1;
  299. for ($j = $degLambda; $j > 0; $j--) {
  300. if ($reg[$j] !== $this->blockSize) {
  301. $reg[$j] = $this->modNn($reg[$j] + $j);
  302. $q ^= $this->alphaTo[$reg[$j]];
  303. }
  304. }
  305. if ($q !== 0) {
  306. // Not a root
  307. continue;
  308. }
  309. // Store root (index-form) and error location number
  310. $root[$count] = $i;
  311. $loc[$count] = $k;
  312. if (++$count === $degLambda) {
  313. break;
  314. }
  315. }
  316. if ($degLambda !== $count) {
  317. // deg(lambda) unequal to number of roots: uncorrectable error detected
  318. return null;
  319. }
  320. // Compute err+eras evaluate poly omega(x) = s(x)*lambda(x) (modulo x**numRoots). In index form. Also find
  321. // deg(omega).
  322. $degOmega = $degLambda - 1;
  323. for ($i = 0; $i <= $degOmega; ++$i) {
  324. $tmp = 0;
  325. for ($j = $i; $j >= 0; --$j) {
  326. if ($syndromes[$i - $j] !== $this->blockSize && $lambda[$j] !== $this->blockSize) {
  327. $tmp ^= $this->alphaTo[$this->modNn($syndromes[$i - $j] + $lambda[$j])];
  328. }
  329. }
  330. $omega[$i] = $this->indexOf[$tmp];
  331. }
  332. // Compute error values in poly-form. num1 = omega(inv(X(l))), num2 = inv(X(l))**(firstRoot-1) and
  333. // den = lambda_pr(inv(X(l))) all in poly form.
  334. for ($j = $count - 1; $j >= 0; --$j) {
  335. $num1 = 0;
  336. for ($i = $degOmega; $i >= 0; $i--) {
  337. if ($omega[$i] !== $this->blockSize) {
  338. $num1 ^= $this->alphaTo[$this->modNn($omega[$i] + $i * $root[$j])];
  339. }
  340. }
  341. $num2 = $this->alphaTo[$this->modNn($root[$j] * ($this->firstRoot - 1) + $this->blockSize)];
  342. $den = 0;
  343. // lambda[i+1] for i even is the formal derivativelambda_pr of lambda[i]
  344. for ($i = min($degLambda, $this->numRoots - 1) & ~1; $i >= 0; $i -= 2) {
  345. if ($lambda[$i + 1] !== $this->blockSize) {
  346. $den ^= $this->alphaTo[$this->modNn($lambda[$i + 1] + $i * $root[$j])];
  347. }
  348. }
  349. // Apply error to data
  350. if ($num1 !== 0 && $loc[$j] >= $this->padding) {
  351. $data[$loc[$j] - $this->padding] = $data[$loc[$j] - $this->padding] ^ (
  352. $this->alphaTo[
  353. $this->modNn(
  354. $this->indexOf[$num1] + $this->indexOf[$num2] + $this->blockSize - $this->indexOf[$den]
  355. )
  356. ]
  357. );
  358. }
  359. }
  360. if (null !== $erasures) {
  361. if (count($erasures) < $count) {
  362. $erasures->setSize($count);
  363. }
  364. for ($i = 0; $i < $count; $i++) {
  365. $erasures[$i] = $loc[$i];
  366. }
  367. }
  368. return $count;
  369. }
  370. /**
  371. * Computes $x % GF_SIZE, where GF_SIZE is 2**GF_BITS - 1, without a slow divide.
  372. */
  373. private function modNn(int $x) : int
  374. {
  375. while ($x >= $this->blockSize) {
  376. $x -= $this->blockSize;
  377. $x = ($x >> $this->symbolSize) + ($x & $this->blockSize);
  378. }
  379. return $x;
  380. }
  381. }