mirror of
https://github.com/claunia/cuetools.net.git
synced 2025-12-16 18:14:25 +00:00
632 lines
28 KiB
C#
632 lines
28 KiB
C#
using System;
|
||
using System.Collections.Generic;
|
||
using System.Text;
|
||
|
||
namespace CUETools.Parity
|
||
{
|
||
public class Galois
|
||
{
|
||
private ushort[] expTbl; // 二重にもつことによりmul, div等を簡略化
|
||
private ushort[] logTbl;
|
||
private int w;
|
||
private int max;
|
||
private int symStart = 0;
|
||
|
||
/**
|
||
* スカラー、ベクターの相互変換テーブルの作成
|
||
*/
|
||
public Galois(int polynomial, int _w)
|
||
{
|
||
w = _w;
|
||
max = (1 << _w) - 1;
|
||
|
||
expTbl = new ushort[max * 2];
|
||
logTbl = new ushort[max + 1];
|
||
|
||
int d = 1;
|
||
for (int i = 0; i < max; i++)
|
||
{
|
||
//if (d == 0)
|
||
// throw new Exception("oops");
|
||
expTbl[i] = expTbl[max + i] = (ushort)d;
|
||
logTbl[d] = (ushort)i;
|
||
d <<= 1;
|
||
if (((d >> _w) & 1) != 0)
|
||
d = (d ^ polynomial) & max;
|
||
}
|
||
}
|
||
|
||
public int Max
|
||
{
|
||
get
|
||
{
|
||
return max;
|
||
}
|
||
}
|
||
|
||
public ushort[] ExpTbl
|
||
{
|
||
get
|
||
{
|
||
return expTbl;
|
||
}
|
||
}
|
||
|
||
public ushort[] LogTbl
|
||
{
|
||
get
|
||
{
|
||
return logTbl;
|
||
}
|
||
}
|
||
|
||
/**
|
||
* スカラー -> ベクター変換
|
||
*
|
||
* @param a int
|
||
* @return int
|
||
*/
|
||
public int toExp(int a)
|
||
{
|
||
return expTbl[a];
|
||
}
|
||
|
||
/**
|
||
* ベクター -> スカラー変換
|
||
*
|
||
* @param a int
|
||
* @return int
|
||
*/
|
||
public int toLog(int a)
|
||
{
|
||
return logTbl[a];
|
||
}
|
||
|
||
/**
|
||
* 誤り位置インデックスの計算
|
||
*
|
||
* @param length int
|
||
* データ長
|
||
* @param a int
|
||
* 誤り位置ベクター
|
||
* @return int
|
||
* 誤り位置インデックス
|
||
*/
|
||
public int toPos(int length, int a)
|
||
{
|
||
return length - 1 - logTbl[a];
|
||
}
|
||
|
||
/**
|
||
* 掛け算
|
||
*
|
||
* @param a int
|
||
* @param b int
|
||
* @return int
|
||
* = a * b
|
||
*/
|
||
public int mul(int a, int b)
|
||
{
|
||
return (a == 0 || b == 0) ? 0 : expTbl[(int)logTbl[a] + logTbl[b]];
|
||
}
|
||
|
||
/**
|
||
* 掛け算
|
||
*
|
||
* @param a int
|
||
* @param b int
|
||
* @return int
|
||
* = a * α^b
|
||
*/
|
||
public int mulExp(int a, int b)
|
||
{
|
||
return (a == 0) ? 0 : expTbl[logTbl[a] + b];
|
||
}
|
||
|
||
/**
|
||
* 割り算
|
||
*
|
||
* @param a int
|
||
* @param b int
|
||
* @return int
|
||
* = a / b
|
||
*/
|
||
public int div(int a, int b)
|
||
{
|
||
return (a == 0) ? 0 : expTbl[logTbl[a] - logTbl[b] + max];
|
||
}
|
||
|
||
/**
|
||
* 割り算
|
||
*
|
||
* @param a int
|
||
* @param b int
|
||
* @return int
|
||
* = a / α^b
|
||
*/
|
||
public int divExp(int a, int b)
|
||
{
|
||
return (a == 0) ? 0 : expTbl[logTbl[a] - b + max];
|
||
}
|
||
|
||
/**
|
||
* 逆数
|
||
*
|
||
* @param a int
|
||
* @return int
|
||
* = 1/a
|
||
*/
|
||
public int inv(int a)
|
||
{
|
||
return expTbl[max - logTbl[a]];
|
||
}
|
||
|
||
public int[] toLog(int[] a)
|
||
{
|
||
var res = new int[a.Length];
|
||
for (int i = 0; i < a.Length; i++)
|
||
res[i] = a[i] == 0 ? - 1 : toLog(a[i]);
|
||
return res;
|
||
}
|
||
|
||
public int[] toLog(ushort[] a)
|
||
{
|
||
var res = new int[a.Length];
|
||
for (int i = 0; i < a.Length; i++)
|
||
res[i] = a[i] == 0 ? -1 : toLog(a[i]);
|
||
return res;
|
||
}
|
||
|
||
public int[] toExp(int[] a)
|
||
{
|
||
var res = new int[a.Length];
|
||
for (int i = 0; i < a.Length; i++)
|
||
res[i] = a[i] == -1 ? 0 : toExp(a[i]);
|
||
return res;
|
||
}
|
||
|
||
public int gfadd(int a, int b)
|
||
{
|
||
var a_exp = a == -1 ? 0 : toExp(a);
|
||
var b_exp = b == -1 ? 0 : toExp(b);
|
||
var res_exp = a_exp ^ b_exp;
|
||
return res_exp == 0 ? -1 : toLog(res_exp);
|
||
}
|
||
|
||
public int[] gfadd(int[] a, int b)
|
||
{
|
||
var res = new int[a.Length];
|
||
var a_exp = toExp(a);
|
||
var b_exp = b == -1 ? 0 : toExp(b);
|
||
for (int i = 0; i < a.Length; i++)
|
||
res[i] = a_exp[i] ^ b_exp;
|
||
return toLog(res);
|
||
}
|
||
|
||
public int[] gfdiff(int[] a)
|
||
{
|
||
//l = length(polynomial);
|
||
//for cc = 2:l
|
||
// %cc-1 represents the power of x
|
||
// if mod(cc-1,2) == 0 %all the even powers are zero because of GF(2)
|
||
// diff(cc-1) = -Inf;
|
||
// else
|
||
// diff(cc-1) = polynomial(cc);
|
||
// end
|
||
//end
|
||
var res = new int[a.Length - 1];
|
||
for (int i = 0; i < res.Length; i++)
|
||
res[i] = (i % 2) == 0 ? a[i + 1] : -1;
|
||
return res;
|
||
}
|
||
|
||
public int gfmul(int a, int b)
|
||
{
|
||
return a < 0 || b < 0 ? -1 : ((a + b) % max);
|
||
}
|
||
|
||
public int gfdiv(int a, int b)
|
||
{
|
||
return a < 0 ? -1 : ((max + a - b) % max);
|
||
}
|
||
|
||
public int gfpow(int value, int p)
|
||
{
|
||
return (value * p) % max;
|
||
}
|
||
|
||
public int gfsubstitute(int[] polynomial, int value, int terms)
|
||
{
|
||
var sum = 0;
|
||
if (value != -1)
|
||
for (int p = 0; p < terms; p++)
|
||
if (polynomial[p] != -1)
|
||
{
|
||
var pow = polynomial[p] + value * p;
|
||
sum ^= expTbl[(pow & max) + (pow >> w)];
|
||
}
|
||
return sum == 0 ? -1 : logTbl[sum];
|
||
}
|
||
|
||
public int[] gfconv(int[] a, int[] b)
|
||
{
|
||
return gfconv(a, b, a.Length + b.Length - 1);
|
||
}
|
||
|
||
public int[] gfconv(int[] a, int[] b, int len)
|
||
{
|
||
var seki = new int[len];
|
||
for (int ia = 0; ia < a.Length; ia++)
|
||
{
|
||
var loga = a[ia];
|
||
if (loga != -1)
|
||
{
|
||
int ib2 = Math.Min(b.Length, len - ia);
|
||
for (int ib = 0; ib < ib2; ib++)
|
||
{
|
||
var logb = b[ib];
|
||
if (logb != -1)
|
||
seki[ia + ib] ^= expTbl[loga + logb]; // = a[ia] * b[ib]
|
||
}
|
||
}
|
||
}
|
||
for (int i = 0; i < len; i++)
|
||
seki[i] = seki[i] == 0 ? -1 : logTbl[seki[i]];
|
||
return seki;
|
||
}
|
||
|
||
public unsafe void gfconv(int* a, int alen, int* b, int blen, int* c, int clen)
|
||
{
|
||
for (int i = 0; i < clen; i++)
|
||
c[i] = 0;
|
||
for (int ia = 0; ia < alen; ia++)
|
||
{
|
||
var loga = a[ia];
|
||
if (loga != -1)
|
||
{
|
||
int ib2 = Math.Min(blen, clen - ia);
|
||
for (int ib = 0; ib < ib2; ib++)
|
||
{
|
||
var logb = b[ib];
|
||
if (logb != -1)
|
||
c[ia + ib] ^= expTbl[loga + logb]; // = a[ia] * b[ib]
|
||
}
|
||
}
|
||
}
|
||
for (int i = 0; i < clen; i++)
|
||
c[i] = c[i] == 0 ? -1 : logTbl[c[i]];
|
||
}
|
||
|
||
public int[] mulPoly(int[] a, int[] b)
|
||
{
|
||
return mulPoly(a, b, a.Length + b.Length - 1);
|
||
}
|
||
|
||
public int[] mulPoly(int[] a, int[] b, int len)
|
||
{
|
||
var res = new int[len];
|
||
mulPoly(res, a, b);
|
||
return res;
|
||
}
|
||
|
||
private static int[] erasures_pos;
|
||
private static int[] erasure_loc_pol;
|
||
private static int[] erasure_diff;
|
||
|
||
public unsafe void Syndrome2Parity(ushort* syndrome, ushort* parity, int npar)
|
||
{
|
||
if (erasure_loc_pol == null)
|
||
{
|
||
var num_erasures = npar;
|
||
|
||
// Compute the erasure locator polynomial:
|
||
erasures_pos = new int[num_erasures];
|
||
for (int x = 0; x < num_erasures; x++)
|
||
erasures_pos[x] = x;
|
||
|
||
//%Compute the erasure-locator polynomial
|
||
// Optimized version
|
||
var erasure_loc_pol_exp = new int[num_erasures + 1];
|
||
erasure_loc_pol_exp[0] = 1;
|
||
for (int i = 0; i < num_erasures; i++)
|
||
for (int x = num_erasures; x > 0; x--)
|
||
erasure_loc_pol_exp[x] ^= mulExp(erasure_loc_pol_exp[x - 1], erasures_pos[i]);
|
||
erasure_loc_pol = toLog(erasure_loc_pol_exp);
|
||
erasure_diff = gfdiff(erasure_loc_pol);
|
||
}
|
||
|
||
// Advance syndrome by npar zeroes (for npar 'erased' parity symbols)
|
||
var S_pol = new int[npar + 1];
|
||
S_pol[0] = -1;
|
||
for (int i = 0; i < npar; i++)
|
||
{
|
||
if (syndrome[i] == 0)
|
||
S_pol[i + 1] = -1;
|
||
else
|
||
{
|
||
var exp = logTbl[syndrome[i]] + npar * i;
|
||
S_pol[i + 1] = (exp & max) + (exp >> w);
|
||
}
|
||
}
|
||
var mod_syn = gfconv(erasure_loc_pol, S_pol, npar + 1);
|
||
|
||
//%Calculate remaining errors (non-erasures)
|
||
//
|
||
//S_M = [];
|
||
//for i = 1:h - num_erasures
|
||
// S_M(i) = mod_syn(i + num_erasures + 1);
|
||
//end
|
||
//flag = 0;
|
||
//if isempty(S_M) == 1
|
||
// flag = 0;
|
||
//else
|
||
// for i = 1:length(S_M)
|
||
// if (S_M(i) ~= -Inf)
|
||
// flag = 1; %Other errors occured in conjunction with erasures
|
||
// end
|
||
// end
|
||
//end
|
||
//%Find error-location polynomial sigma (Berlekamp's iterative algorithm -
|
||
//if (flag == 1)
|
||
//{
|
||
// ...
|
||
//}
|
||
//else
|
||
//{
|
||
// sigma = 0;
|
||
// comp_error_locs = [];
|
||
//}
|
||
|
||
#if kjsljdf
|
||
var sigma = new int[1] { 0 };
|
||
var omega = gfconv(sigma, gfadd(mod_syn, 0), npar + 1);
|
||
var tsi = gfconv(sigma, erasure_loc_pol);
|
||
var tsi_diff = gfdiff(tsi);
|
||
//var e_e_places = [erasures_pos comp_error_locs];
|
||
#else
|
||
var omega = mod_syn;
|
||
var tsi_diff = erasure_diff;
|
||
var e_e_places = erasures_pos;
|
||
#endif
|
||
|
||
//%Calculate the error magnitudes
|
||
//%Substitute the inverse into sigma_diff
|
||
//var ERR = new int[e_e_places.Length];
|
||
for (int ii = 0; ii < e_e_places.Length; ii++)
|
||
{
|
||
var point = max - e_e_places[ii];
|
||
var ERR_DEN = gfsubstitute(tsi_diff, point, tsi_diff.Length);
|
||
var ERR_NUM = gfsubstitute(omega, point, omega.Length);
|
||
// Additional +ii because we use slightly different syndrome
|
||
var pow = ERR_NUM + e_e_places[ii] + ii + max - ERR_DEN;
|
||
|
||
//ERR[ii] = ERR_NUM == -1 ? 0 : expTbl[(pow & max) + (pow >> w)];
|
||
parity[npar - 1 - ii] = ERR_NUM == -1 ? (ushort)0 : expTbl[(pow & max) + (pow >> w)];
|
||
}
|
||
}
|
||
|
||
/**
|
||
* 数式の掛け算
|
||
*
|
||
* @param seki int[]
|
||
* seki = a * b
|
||
* @param a int[]
|
||
* @param b int[]
|
||
*/
|
||
public void mulPoly(int[] seki, int[] a, int[] b)
|
||
{
|
||
Array.Clear(seki, 0, seki.Length);
|
||
for (int ia = 0; ia < a.Length; ia++)
|
||
{
|
||
if (a[ia] != 0)
|
||
{
|
||
int loga = logTbl[a[ia]];
|
||
int ib2 = Math.Min(b.Length, seki.Length - ia);
|
||
for (int ib = 0; ib < ib2; ib++)
|
||
{
|
||
if (b[ib] != 0)
|
||
{
|
||
seki[ia + ib] ^= expTbl[loga + logTbl[b[ib]]]; // = a[ia] * b[ib]
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
public unsafe void mulPoly(int* seki, int* a, int* b, int lenS, int lenA, int lenB)
|
||
{
|
||
for (int i = 0; i < lenS; i++)
|
||
seki[i] = 0;
|
||
for (int ia = 0; ia < lenA; ia++)
|
||
{
|
||
if (a[ia] != 0)
|
||
{
|
||
int loga = logTbl[a[ia]];
|
||
int ib2 = Math.Min(lenB, lenS - ia);
|
||
for (int ib = 0; ib < ib2; ib++)
|
||
{
|
||
if (b[ib] != 0)
|
||
{
|
||
seki[ia + ib] ^= expTbl[loga + logTbl[b[ib]]]; // = a[ia] * b[ib]
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
/**
|
||
* 生成多項式配列の作成
|
||
* G(x)=Π[k=0,n-1](x + α^k)
|
||
* encodeGxの添え字と次数の並びが逆なのに注意
|
||
* encodeGx[0] = x^(npar - 1)の項
|
||
* encodeGx[1] = x^(npar - 2)の項
|
||
* ...
|
||
* encodeGx[npar - 1] = x^0の項
|
||
*/
|
||
public int[] makeEncodeGx(int npar)
|
||
{
|
||
int[] encodeGx = new int[npar];
|
||
encodeGx[npar - 1] = 1;
|
||
for (int i = 0, kou = symStart; i < npar; i++, kou++)
|
||
{
|
||
int ex = toExp(kou); // ex = α^kou
|
||
// (x + α^kou)を掛る
|
||
for (int j = 0; j < npar - 1; j++)
|
||
{
|
||
// 現在の項 * α^kou + 一つ下の次数の項
|
||
encodeGx[j] = mul(encodeGx[j], ex) ^ encodeGx[j + 1];
|
||
}
|
||
encodeGx[npar - 1] = mul(encodeGx[npar - 1], ex);// 最下位項の計算
|
||
}
|
||
return encodeGx;
|
||
}
|
||
|
||
public int[] makeEncodeGxLog(int npar)
|
||
{
|
||
int[] encodeGx = makeEncodeGx(npar);
|
||
for (int i = 0; i < npar; i++)
|
||
{
|
||
if (encodeGx[i] == 0)
|
||
throw new Exception("0 in encodeGx");
|
||
encodeGx[i] = toLog(encodeGx[i]);
|
||
}
|
||
return encodeGx;
|
||
}
|
||
|
||
/// <summary>
|
||
/// parityTable[xx, 0, i] = mul(00xx, encodeGx[i])
|
||
/// parityTable[xx, 1, i] = mul(xx00, encodeGx[i])
|
||
/// </summary>
|
||
/// <param name="npar"></param>
|
||
/// <returns></returns>
|
||
public ushort[,,] makeEncodeTable(int npar)
|
||
{
|
||
var loggx = makeEncodeGxLog(npar);
|
||
var parityTable = new ushort[256, 2, npar];
|
||
for (int i = 0; i < npar; i++)
|
||
{
|
||
parityTable[0, 0, i] = 0;
|
||
parityTable[0, 1, i] = 0;
|
||
}
|
||
for (int ib = 1; ib < 256; ib++)
|
||
{
|
||
int logib0 = LogTbl[ib];
|
||
int logib1 = LogTbl[ib << 8];
|
||
for (int i = 0; i < npar; i++)
|
||
{
|
||
parityTable[ib, 0, i] = ExpTbl[logib0 + loggx[i]];
|
||
parityTable[ib, 1, i] = ExpTbl[logib1 + loggx[i]];
|
||
}
|
||
}
|
||
return parityTable;
|
||
}
|
||
|
||
/// <summary>
|
||
/// parityTable[xx, 0, i] = mul(00xx, α^i)
|
||
/// parityTable[xx, 1, i] = mul(xx00, α^i)
|
||
/// </summary>
|
||
/// <param name="npar"></param>
|
||
/// <returns></returns>
|
||
public ushort[,,] makeDecodeTable(int npar)
|
||
{
|
||
var parityTable = new ushort[256, 2, npar];
|
||
for (int i = 0; i < npar; i++)
|
||
{
|
||
parityTable[0, 0, i] = 0;
|
||
parityTable[0, 1, i] = 0;
|
||
}
|
||
for (int ib = 1; ib < 256; ib++)
|
||
{
|
||
int logib0 = LogTbl[ib];
|
||
int logib1 = LogTbl[ib << 8];
|
||
for (int i = 0; i < npar; i++)
|
||
{
|
||
parityTable[ib, 0, i] = ExpTbl[logib0 + i];
|
||
parityTable[ib, 1, i] = ExpTbl[logib1 + i];
|
||
}
|
||
}
|
||
return parityTable;
|
||
}
|
||
|
||
/**
|
||
* シンドロームの計算
|
||
* @param data int[]
|
||
* 入力データ配列
|
||
* @param length int
|
||
* データ長
|
||
* @param syn int[]
|
||
* (x - α^0) (x - α^1) (x - α^2) ...のシンドローム
|
||
* @return boolean
|
||
* true: シンドロームは総て0
|
||
*/
|
||
public bool calcSyndrome(byte[] data, int length, int[] syn)
|
||
{
|
||
int hasErr = 0;
|
||
for (int i = 0; i < syn.Length; i++)
|
||
{
|
||
int wk = 0;
|
||
for (int idx = 0; idx < length; idx++)
|
||
{
|
||
//wk = data[idx] ^ ((wk == 0) ? 0 : expTbl[logTbl[wk] + i + symStart]); // wk = data + wk * α^i
|
||
wk = data[idx] ^ ((wk == 0) ? 0 : expTbl[logTbl[wk] + i]); // wk = data + wk * α^i
|
||
}
|
||
syn[i] = wk;
|
||
hasErr |= wk;
|
||
}
|
||
return hasErr == 0;
|
||
}
|
||
|
||
/**
|
||
* シンドロームの計算
|
||
* @param data int[]
|
||
* 入力データ配列
|
||
* @param length int
|
||
* データ長
|
||
* @param syn int[]
|
||
* (x - α^0) (x - α^1) (x - α^2) ...のシンドローム
|
||
* @return boolean
|
||
* true: シンドロームは総て0
|
||
*/
|
||
public unsafe bool calcSyndrome(ushort* data, int length, int[] syn)
|
||
{
|
||
int hasErr = 0;
|
||
for (int i = 0; i < syn.Length; i++)
|
||
{
|
||
int wk = 0;
|
||
for (int idx = 0; idx < length; idx++)
|
||
{
|
||
//wk = data[idx] ^ ((wk == 0) ? 0 : expTbl[logTbl[wk] + i + symStart]); // wk = data + wk * α^i
|
||
wk = data[idx] ^ ((wk == 0) ? 0 : expTbl[logTbl[wk] + i]); // wk = data + wk * α^i
|
||
}
|
||
syn[i] = wk;
|
||
hasErr |= wk;
|
||
}
|
||
return hasErr == 0;
|
||
}
|
||
}
|
||
|
||
public class Galois81D: Galois
|
||
{
|
||
public const int POLYNOMIAL = 0x1d;
|
||
|
||
public static Galois81D instance = new Galois81D();
|
||
|
||
Galois81D()
|
||
: base(POLYNOMIAL, 8)
|
||
{
|
||
}
|
||
}
|
||
|
||
public class Galois16 : Galois
|
||
{
|
||
public const int POLYNOMIAL = 0x1100B;
|
||
|
||
public static Galois16 instance = new Galois16();
|
||
|
||
Galois16()
|
||
: base(POLYNOMIAL, 16)
|
||
{
|
||
}
|
||
}
|
||
}
|