Skip to content

Instantly share code, notes, and snippets.

@lifthrasiir

lifthrasiir/rs_ccsds.py

Created November 14, 2023 18:21
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save lifthrasiir/6b8900b0d55bba8f1745fdb82cf1d60d to your computer and use it in GitHub Desktop.
Save lifthrasiir/6b8900b0d55bba8f1745fdb82cf1d60d to your computer and use it in GitHub Desktop.
Download ZIP
Optimized Python version of crozone/ReedSolomonCCSDS
Raw
rs_ccsds.py
# derived from https://github.com/crozone/ReedSolomonCCSDS
# should be more performant alternative to https://github.com/kplabs-pl/reed-solomon-ccsds
from dataclasses import dataclass
from typing import Optional, List
ALPHA_TO = [
0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80,0x87,0x89,0x95,0xad,0xdd,0x3d,0x7a,0xf4,
0x6f,0xde,0x3b,0x76,0xec,0x5f,0xbe,0xfb,0x71,0xe2,0x43,0x86,0x8b,0x91,0xa5,0xcd,
0x1d,0x3a,0x74,0xe8,0x57,0xae,0xdb,0x31,0x62,0xc4,0x0f,0x1e,0x3c,0x78,0xf0,0x67,
0xce,0x1b,0x36,0x6c,0xd8,0x37,0x6e,0xdc,0x3f,0x7e,0xfc,0x7f,0xfe,0x7b,0xf6,0x6b,
0xd6,0x2b,0x56,0xac,0xdf,0x39,0x72,0xe4,0x4f,0x9e,0xbb,0xf1,0x65,0xca,0x13,0x26,
0x4c,0x98,0xb7,0xe9,0x55,0xaa,0xd3,0x21,0x42,0x84,0x8f,0x99,0xb5,0xed,0x5d,0xba,
0xf3,0x61,0xc2,0x03,0x06,0x0c,0x18,0x30,0x60,0xc0,0x07,0x0e,0x1c,0x38,0x70,0xe0,
0x47,0x8e,0x9b,0xb1,0xe5,0x4d,0x9a,0xb3,0xe1,0x45,0x8a,0x93,0xa1,0xc5,0x0d,0x1a,
0x34,0x68,0xd0,0x27,0x4e,0x9c,0xbf,0xf9,0x75,0xea,0x53,0xa6,0xcb,0x11,0x22,0x44,
0x88,0x97,0xa9,0xd5,0x2d,0x5a,0xb4,0xef,0x59,0xb2,0xe3,0x41,0x82,0x83,0x81,0x85,
0x8d,0x9d,0xbd,0xfd,0x7d,0xfa,0x73,0xe6,0x4b,0x96,0xab,0xd1,0x25,0x4a,0x94,0xaf,
0xd9,0x35,0x6a,0xd4,0x2f,0x5e,0xbc,0xff,0x79,0xf2,0x63,0xc6,0x0b,0x16,0x2c,0x58,
0xb0,0xe7,0x49,0x92,0xa3,0xc1,0x05,0x0a,0x14,0x28,0x50,0xa0,0xc7,0x09,0x12,0x24,
0x48,0x90,0xa7,0xc9,0x15,0x2a,0x54,0xa8,0xd7,0x29,0x52,0xa4,0xcf,0x19,0x32,0x64,
0xc8,0x17,0x2e,0x5c,0xb8,0xf7,0x69,0xd2,0x23,0x46,0x8c,0x9f,0xb9,0xf5,0x6d,0xda,
0x33,0x66,0xcc,0x1f,0x3e,0x7c,0xf8,0x77,0xee,0x5b,0xb6,0xeb,0x51,0xa2,0xc3,0x00
]
INDEX_OF = [
0xff,0x00,0x01,0x63,0x02,0xc6,0x64,0x6a,0x03,0xcd,0xc7,0xbc,0x65,0x7e,0x6b,0x2a,
0x04,0x8d,0xce,0x4e,0xc8,0xd4,0xbd,0xe1,0x66,0xdd,0x7f,0x31,0x6c,0x20,0x2b,0xf3,
0x05,0x57,0x8e,0xe8,0xcf,0xac,0x4f,0x83,0xc9,0xd9,0xd5,0x41,0xbe,0x94,0xe2,0xb4,
0x67,0x27,0xde,0xf0,0x80,0xb1,0x32,0x35,0x6d,0x45,0x21,0x12,0x2c,0x0d,0xf4,0x38,
0x06,0x9b,0x58,0x1a,0x8f,0x79,0xe9,0x70,0xd0,0xc2,0xad,0xa8,0x50,0x75,0x84,0x48,
0xca,0xfc,0xda,0x8a,0xd6,0x54,0x42,0x24,0xbf,0x98,0x95,0xf9,0xe3,0x5e,0xb5,0x15,
0x68,0x61,0x28,0xba,0xdf,0x4c,0xf1,0x2f,0x81,0xe6,0xb2,0x3f,0x33,0xee,0x36,0x10,
0x6e,0x18,0x46,0xa6,0x22,0x88,0x13,0xf7,0x2d,0xb8,0x0e,0x3d,0xf5,0xa4,0x39,0x3b,
0x07,0x9e,0x9c,0x9d,0x59,0x9f,0x1b,0x08,0x90,0x09,0x7a,0x1c,0xea,0xa0,0x71,0x5a,
0xd1,0x1d,0xc3,0x7b,0xae,0x0a,0xa9,0x91,0x51,0x5b,0x76,0x72,0x85,0xa1,0x49,0xeb,
0xcb,0x7c,0xfd,0xc4,0xdb,0x1e,0x8b,0xd2,0xd7,0x92,0x55,0xaa,0x43,0x0b,0x25,0xaf,
0xc0,0x73,0x99,0x77,0x96,0x5c,0xfa,0x52,0xe4,0xec,0x5f,0x4a,0xb6,0xa2,0x16,0x86,
0x69,0xc5,0x62,0xfe,0x29,0x7d,0xbb,0xcc,0xe0,0xd3,0x4d,0x8c,0xf2,0x1f,0x30,0xdc,
0x82,0xab,0xe7,0x56,0xb3,0x93,0x40,0xd8,0x34,0xb0,0xef,0x26,0x37,0x0c,0x11,0x44,
0x6f,0x78,0x19,0x9a,0x47,0x74,0xa7,0xc1,0x23,0x53,0x89,0xfb,0x14,0x5d,0xf8,0x97,
0x2e,0x4b,0xb9,0x60,0x0f,0xed,0x3e,0xe5,0xf6,0x87,0xa5,0x17,0x3a,0xa3,0x3c,0xb7
]
TAL_TO_DUAL_BASIS = [
0x00,0x7b,0xaf,0xd4,0x99,0xe2,0x36,0x4d,0xfa,0x81,0x55,0x2e,0x63,0x18,0xcc,0xb7,
0x86,0xfd,0x29,0x52,0x1f,0x64,0xb0,0xcb,0x7c,0x07,0xd3,0xa8,0xe5,0x9e,0x4a,0x31,
0xec,0x97,0x43,0x38,0x75,0x0e,0xda,0xa1,0x16,0x6d,0xb9,0xc2,0x8f,0xf4,0x20,0x5b,
0x6a,0x11,0xc5,0xbe,0xf3,0x88,0x5c,0x27,0x90,0xeb,0x3f,0x44,0x09,0x72,0xa6,0xdd,
0xef,0x94,0x40,0x3b,0x76,0x0d,0xd9,0xa2,0x15,0x6e,0xba,0xc1,0x8c,0xf7,0x23,0x58,
0x69,0x12,0xc6,0xbd,0xf0,0x8b,0x5f,0x24,0x93,0xe8,0x3c,0x47,0x0a,0x71,0xa5,0xde,
0x03,0x78,0xac,0xd7,0x9a,0xe1,0x35,0x4e,0xf9,0x82,0x56,0x2d,0x60,0x1b,0xcf,0xb4,
0x85,0xfe,0x2a,0x51,0x1c,0x67,0xb3,0xc8,0x7f,0x04,0xd0,0xab,0xe6,0x9d,0x49,0x32,
0x8d,0xf6,0x22,0x59,0x14,0x6f,0xbb,0xc0,0x77,0x0c,0xd8,0xa3,0xee,0x95,0x41,0x3a,
0x0b,0x70,0xa4,0xdf,0x92,0xe9,0x3d,0x46,0xf1,0x8a,0x5e,0x25,0x68,0x13,0xc7,0xbc,
0x61,0x1a,0xce,0xb5,0xf8,0x83,0x57,0x2c,0x9b,0xe0,0x34,0x4f,0x02,0x79,0xad,0xd6,
0xe7,0x9c,0x48,0x33,0x7e,0x05,0xd1,0xaa,0x1d,0x66,0xb2,0xc9,0x84,0xff,0x2b,0x50,
0x62,0x19,0xcd,0xb6,0xfb,0x80,0x54,0x2f,0x98,0xe3,0x37,0x4c,0x01,0x7a,0xae,0xd5,
0xe4,0x9f,0x4b,0x30,0x7d,0x06,0xd2,0xa9,0x1e,0x65,0xb1,0xca,0x87,0xfc,0x28,0x53,
0x8e,0xf5,0x21,0x5a,0x17,0x6c,0xb8,0xc3,0x74,0x0f,0xdb,0xa0,0xed,0x96,0x42,0x39,
0x08,0x73,0xa7,0xdc,0x91,0xea,0x3e,0x45,0xf2,0x89,0x5d,0x26,0x6b,0x10,0xc4,0xbf
]
TAL_TO_CONVENTIONAL = [
0x00,0xcc,0xac,0x60,0x79,0xb5,0xd5,0x19,0xf0,0x3c,0x5c,0x90,0x89,0x45,0x25,0xe9,
0xfd,0x31,0x51,0x9d,0x84,0x48,0x28,0xe4,0x0d,0xc1,0xa1,0x6d,0x74,0xb8,0xd8,0x14,
0x2e,0xe2,0x82,0x4e,0x57,0x9b,0xfb,0x37,0xde,0x12,0x72,0xbe,0xa7,0x6b,0x0b,0xc7,
0xd3,0x1f,0x7f,0xb3,0xaa,0x66,0x06,0xca,0x23,0xef,0x8f,0x43,0x5a,0x96,0xf6,0x3a,
0x42,0x8e,0xee,0x22,0x3b,0xf7,0x97,0x5b,0xb2,0x7e,0x1e,0xd2,0xcb,0x07,0x67,0xab,
0xbf,0x73,0x13,0xdf,0xc6,0x0a,0x6a,0xa6,0x4f,0x83,0xe3,0x2f,0x36,0xfa,0x9a,0x56,
0x6c,0xa0,0xc0,0x0c,0x15,0xd9,0xb9,0x75,0x9c,0x50,0x30,0xfc,0xe5,0x29,0x49,0x85,
0x91,0x5d,0x3d,0xf1,0xe8,0x24,0x44,0x88,0x61,0xad,0xcd,0x01,0x18,0xd4,0xb4,0x78,
0xc5,0x09,0x69,0xa5,0xbc,0x70,0x10,0xdc,0x35,0xf9,0x99,0x55,0x4c,0x80,0xe0,0x2c,
0x38,0xf4,0x94,0x58,0x41,0x8d,0xed,0x21,0xc8,0x04,0x64,0xa8,0xb1,0x7d,0x1d,0xd1,
0xeb,0x27,0x47,0x8b,0x92,0x5e,0x3e,0xf2,0x1b,0xd7,0xb7,0x7b,0x62,0xae,0xce,0x02,
0x16,0xda,0xba,0x76,0x6f,0xa3,0xc3,0x0f,0xe6,0x2a,0x4a,0x86,0x9f,0x53,0x33,0xff,
0x87,0x4b,0x2b,0xe7,0xfe,0x32,0x52,0x9e,0x77,0xbb,0xdb,0x17,0x0e,0xc2,0xa2,0x6e,
0x7a,0xb6,0xd6,0x1a,0x03,0xcf,0xaf,0x63,0x8a,0x46,0x26,0xea,0xf3,0x3f,0x5f,0x93,
0xa9,0x65,0x05,0xc9,0xd0,0x1c,0x7c,0xb0,0x59,0x95,0xf5,0x39,0x20,0xec,0x8c,0x40,
0x54,0x98,0xf8,0x34,0x2d,0xe1,0x81,0x4d,0xa4,0x68,0x08,0xc4,0xdd,0x11,0x71,0xbd
]
BLOCK_LENGTH = 255 # == N
DATA_LENGTH = 255 - 32 # == N - R
PARITY_LENGTH = 32 # == R
# 1. if 0 <= x, y < N, then ALPHA_TO_EX[x + y] == ALPHA_TO[(x + y) % N]
# 2. if 0 <= x <= N and 0 <= y < N, then
# ALPHA_TO_EX[INDEX_OF_EX[x] + y] == (ALPHA_TO[(INDEX_OF[x] + y) % N] if x != N else 0)
ALPHA_TO_EX = ALPHA_TO[:255] * 2 + [0] * (2 * 255)
INDEX_OF_EX = [2 * 255] + INDEX_OF[1:]
@dataclass
class Encoded:
data: List[int]
parity: List[int]
def encode_block(data: List[int], dual_basis: bool = False) -> Encoded:
assert len(data) == DATA_LENGTH
if dual_basis:
data = [TAL_TO_CONVENTIONAL[i] for i in data]
p0 = p1 = p2 = p3 = p4 = p5 = p6 = p7 = p8 = p9 = \
p10 = p11 = p12 = p13 = p14 = p15 = p16 = p17 = p18 = p19 = \
p20 = p21 = p22 = p23 = p24 = p25 = p26 = p27 = p28 = p29 = p30 = p31 = 0
for b in data:
feedback = INDEX_OF[b ^ p0]
if feedback != 255:
p0 = p1 ^ ALPHA_TO_EX[feedback + 0xf9] # GENPOLY[31]
p1 = p2 ^ ALPHA_TO_EX[feedback + 0x3b] # GENPOLY[30]
p2 = p3 ^ ALPHA_TO_EX[feedback + 0x42] # ...
p3 = p4 ^ ALPHA_TO_EX[feedback + 0x04]
p4 = p5 ^ ALPHA_TO_EX[feedback + 0x2b]
p5 = p6 ^ ALPHA_TO_EX[feedback + 0x7e]
p6 = p7 ^ ALPHA_TO_EX[feedback + 0xfb]
p7 = p8 ^ ALPHA_TO_EX[feedback + 0x61]
p8 = p9 ^ ALPHA_TO_EX[feedback + 0x1e]
p9 = p10 ^ ALPHA_TO_EX[feedback + 0x03]
p10 = p11 ^ ALPHA_TO_EX[feedback + 0xd5]
p11 = p12 ^ ALPHA_TO_EX[feedback + 0x32]
p12 = p13 ^ ALPHA_TO_EX[feedback + 0x42]
p13 = p14 ^ ALPHA_TO_EX[feedback + 0xaa]
p14 = p15 ^ ALPHA_TO_EX[feedback + 0x05]
p15 = p16 ^ ALPHA_TO_EX[feedback + 0x18]
p16 = p17 ^ ALPHA_TO_EX[feedback + 0x05]
p17 = p18 ^ ALPHA_TO_EX[feedback + 0xaa]
p18 = p19 ^ ALPHA_TO_EX[feedback + 0x42]
p19 = p20 ^ ALPHA_TO_EX[feedback + 0x32]
p20 = p21 ^ ALPHA_TO_EX[feedback + 0xd5]
p21 = p22 ^ ALPHA_TO_EX[feedback + 0x03]
p22 = p23 ^ ALPHA_TO_EX[feedback + 0x1e]
p23 = p24 ^ ALPHA_TO_EX[feedback + 0x61]
p24 = p25 ^ ALPHA_TO_EX[feedback + 0xfb]
p25 = p26 ^ ALPHA_TO_EX[feedback + 0x7e]
p26 = p27 ^ ALPHA_TO_EX[feedback + 0x2b]
p27 = p28 ^ ALPHA_TO_EX[feedback + 0x04]
p28 = p29 ^ ALPHA_TO_EX[feedback + 0x42]
p29 = p30 ^ ALPHA_TO_EX[feedback + 0x3b] # ...
p30 = p31 ^ ALPHA_TO_EX[feedback + 0xf9] # GENPOLY[1]
p31 = ALPHA_TO_EX[feedback] # GENPOLY[0] == 0
else:
p0 = p1; p1 = p2; p2 = p3; p3 = p4; p4 = p5
p5 = p6; p6 = p7; p7 = p8; p8 = p9; p9 = p10
p10 = p11; p11 = p12; p12 = p13; p13 = p14; p14 = p15
p15 = p16; p16 = p17; p17 = p18; p18 = p19; p19 = p20
p20 = p21; p21 = p22; p22 = p23; p23 = p24; p24 = p25
p25 = p26; p26 = p27; p27 = p28; p28 = p29; p29 = p30
p30 = p31; p31 = 0
parity = [
p0, p1, p2, p3, p4, p5, p6, p7, p8, p9,
p10, p11, p12, p13, p14, p15, p16, p17, p18, p19,
p20, p21, p22, p23, p24, p25, p26, p27, p28, p29, p30, p31,
]
if dual_basis:
data = [TAL_TO_DUAL_BASIS[i] for i in data]
parity = [TAL_TO_DUAL_BASIS[i] for i in parity]
return Encoded(data=data, parity=parity)
def decode_block(
block: Encoded,
erasures: Optional[List[int]] = None,
dual_basis: bool = False,
) -> List[int]:
assert len(block.data) == DATA_LENGTH
assert len(block.parity) == PARITY_LENGTH
erasures = erasures or ()
block = block.data + block.parity
if dual_basis:
block = [TAL_TO_CONVENTIONAL[i] for i in block]
s0 = s1 = s2 = s3 = s4 = s5 = s6 = s7 = s8 = s9 = \
s10 = s11 = s12 = s13 = s14 = s15 = s16 = s17 = s18 = s19 = \
s20 = s21 = s22 = s23 = s24 = s25 = s26 = s27 = s28 = s29 = s30 = s31 = block[0]
for b in block[1:]:
s0 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s0] + (112 + 0) * 11 % 255]
s1 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s1] + (112 + 1) * 11 % 255]
s2 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s2] + (112 + 2) * 11 % 255]
s3 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s3] + (112 + 3) * 11 % 255]
s4 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s4] + (112 + 4) * 11 % 255]
s5 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s5] + (112 + 5) * 11 % 255]
s6 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s6] + (112 + 6) * 11 % 255]
s7 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s7] + (112 + 7) * 11 % 255]
s8 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s8] + (112 + 8) * 11 % 255]
s9 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s9] + (112 + 9) * 11 % 255]
s10 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s10] + (112 + 10) * 11 % 255]
s11 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s11] + (112 + 11) * 11 % 255]
s12 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s12] + (112 + 12) * 11 % 255]
s13 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s13] + (112 + 13) * 11 % 255]
s14 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s14] + (112 + 14) * 11 % 255]
s15 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s15] + (112 + 15) * 11 % 255]
s16 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s16] + (112 + 16) * 11 % 255]
s17 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s17] + (112 + 17) * 11 % 255]
s18 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s18] + (112 + 18) * 11 % 255]
s19 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s19] + (112 + 19) * 11 % 255]
s20 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s20] + (112 + 20) * 11 % 255]
s21 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s21] + (112 + 21) * 11 % 255]
s22 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s22] + (112 + 22) * 11 % 255]
s23 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s23] + (112 + 23) * 11 % 255]
s24 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s24] + (112 + 24) * 11 % 255]
s25 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s25] + (112 + 25) * 11 % 255]
s26 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s26] + (112 + 26) * 11 % 255]
s27 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s27] + (112 + 27) * 11 % 255]
s28 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s28] + (112 + 28) * 11 % 255]
s29 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s29] + (112 + 29) * 11 % 255]
s30 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s30] + (112 + 30) * 11 % 255]
s31 = b ^ ALPHA_TO_EX[INDEX_OF_EX[s31] + (112 + 31) * 11 % 255]
syndromes = [
s0, s1, s2, s3, s4, s5, s6, s7, s8, s9,
s10, s11, s12, s13, s14, s15, s16, s17, s18, s19,
s20, s21, s22, s23, s24, s25, s26, s27, s28, s29, s30, s31,
]
if any(syndromes):
N = 255
R = 32
syndromes = [INDEX_OF[s] for s in syndromes]
lam = [1] + [0] * R
if erasures:
lam[1] = ALPHA_TO[(11 * (N - 1 - erasures[0])) % N]
for i in range(1, len(erasures)):
u = (11 * (N - 1 - erasures[i])) % N
for j in range(i + 1, 0, -1):
lam[j] ^= ALPHA_TO_EX[u + INDEX_OF_EX[lam[j - 1]]]
b = [INDEX_OF[l] for l in lam]
r = len(erasures) + 1
el = len(erasures)
while r <= R:
discr = 0
for i in range(r):
if syndromes[r - i - 1] != N:
discr ^= ALPHA_TO_EX[INDEX_OF_EX[lam[i]] + syndromes[r - i - 1]]
discr = INDEX_OF[discr]
if discr == N:
# B(x) = x * B(x)
b = [N] + b[:-1]
else:
# T(x) = lambda(x) - discr_r * x * B(x)
t = lam[:]
for i in range(R):
if b[i] != N: t[i + 1] ^= ALPHA_TO_EX[discr + b[i]]
if 2 * el <= r + len(erasures) - 1:
el = r + len(erasures) - el
# B(x) = discr_r^(-1) * lambda(x)
b = [N if l == 0 else (INDEX_OF[l] - discr) % N for l in lam]
else:
# B(x) = x * B(x)
b = [N] + b[:-1]
lam = t
r += 1
lam = [INDEX_OF[l] for l in lam]
deg_lam = 0
for i in range(R + 1):
if lam[i] != N: deg_lam = i
reg = lam[:]
roots = []
for i in range(1, N + 1):
q = 1
for j in range(1, deg_lam + 1): # can skip reg[0] == lam[0] == 0
if reg[j] != N:
reg[j] = (reg[j] + j) % N
q ^= ALPHA_TO[reg[j]]
if q == 0:
roots.append(i)
if len(roots) == deg_lam: break
else:
if len(roots) != deg_lam:
raise RuntimeError('uncorrectable error')
omega = []
deg_omega = 0
for i in range(R):
tmp = 0
for j in range(min(deg_lam, i), -1, -1):
if syndromes[i - j] != N and lam[j] != N:
tmp ^= ALPHA_TO_EX[syndromes[i - j] + lam[j]]
if tmp != 0: deg_omega = i
omega.append(INDEX_OF[tmp])
omega.append(N)
while roots:
root = roots.pop()
num1 = 0
for i in range(deg_omega, -1, -1):
if omega[i] != N:
num1 ^= ALPHA_TO[(omega[i] + i * root) % N]
num2 = ALPHA_TO[root * (112 - 1) % N]
den = 0
start = min(deg_lam, R - 1)
for i in reversed(range(0, start, 2)):
if lam[i + 1] != N:
den ^= ALPHA_TO[(lam[i + 1] + i * root) % N]
if den == 0:
raise RuntimeError('uncorrectable error')
if num1 != 0:
loc = (root * 116 - 1) % N
block[loc] ^= ALPHA_TO[(INDEX_OF[num1] + INDEX_OF[num2] - INDEX_OF[den]) % N]
#print('located:', sorted(loc)) # TODO API
if dual_basis:
block = [TAL_TO_DUAL_BASIS[i] for i in block]
return block[:DATA_LENGTH]
def _self_test():
DATA = [
0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08,0x09,0x0a,0x0b,0x0c,0x0d,0x0e,0x0f,
0x10,0x11,0x12,0x13,0x14,0x15,0x16,0x17,0x18,0x19,0x1a,0x1b,0x1c,0x1d,0x1e,0x1f,
0x20,0x21,0x22,0x23,0x24,0x25,0x26,0x27,0x28,0x29,0x2a,0x2b,0x2c,0x2d,0x2e,0x2f,
0x30,0x31,0x32,0x33,0x34,0x35,0x36,0x37,0x38,0x39,0x3a,0x3b,0x3c,0x3d,0x3e,0x3f,
0x40,0x41,0x42,0x43,0x44,0x45,0x46,0x47,0x48,0x49,0x4a,0x4b,0x4c,0x4d,0x4e,0x4f,
0x50,0x51,0x52,0x53,0x54,0x55,0x56,0x57,0x58,0x59,0x5a,0x5b,0x5c,0x5d,0x5e,0x5f,
0x60,0x61,0x62,0x63,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x6b,0x6c,0x6d,0x6e,0x6f,
0x70,0x71,0x72,0x73,0x74,0x75,0x76,0x77,0x78,0x79,0x7a,0x7b,0x7c,0x7d,0x7e,0x7f,
0x80,0x81,0x82,0x83,0x84,0x85,0x86,0x87,0x88,0x89,0x8a,0x8b,0x8c,0x8d,0x8e,0x8f,
0x90,0x91,0x92,0x93,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0x9b,0x9c,0x9d,0x9e,0x9f,
0xa0,0xa1,0xa2,0xa3,0xa4,0xa5,0xa6,0xa7,0xa8,0xa9,0xaa,0xab,0xac,0xad,0xae,0xaf,
0xb0,0xb1,0xb2,0xb3,0xb4,0xb5,0xb6,0xb7,0xb8,0xb9,0xba,0xbb,0xbc,0xbd,0xbe,0xbf,
0xc0,0xc1,0xc2,0xc3,0xc4,0xc5,0xc6,0xc7,0xc8,0xc9,0xca,0xcb,0xcc,0xcd,0xce,0xcf,
0xd0,0xd1,0xd2,0xd3,0xd4,0xd5,0xd6,0xd7,0xd8,0xd9,0xda,0xdb,0xdc,0xdd,0xde
]
PARITY_CONVENTIONAL = [
0x2f,0xbd,0x4f,0xb4,0x74,0x84,0x94,0xb9,0xac,0xd5,0x54,0x62,0x72,0x12,0xee,0xb3,
0xeb,0xed,0x41,0x19,0x1d,0xe1,0xd3,0x63,0x20,0xea,0x49,0x29,0x0b,0x25,0xab,0xcf,
]
PARITY_DUAL_BASIS = [
0x4F,0xFB,0x92,0xDD,0x55,0x7E,0xC6,0x7F,0x27,0xFB,0x89,0x82,0xCF,0x58,0xF8,0xFD,
0x02,0x8A,0xD1,0x17,0xFC,0xEF,0x6B,0x27,0x93,0xD0,0x41,0x88,0x26,0x57,0x86,0x51,
]
DATA_ERROR_MASK = [
0x58,0x00,0xA3,0x00,0x00,0xCD,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x0D,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0xCA,0x96,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x1B,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0xA2,0x00,0xAC,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0xB9,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0xE5,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x94,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0xC3,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x97,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00
]
PARITY_ERROR_MASK = [
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x7A,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x29,0x00,0x00,0x00,0x00,
]
#print('expected:', [i for i, x in enumerate(DATA_ERROR_MASK + PARITY_ERROR_MASK) if x])
encoded = encode_block(DATA)
assert encoded.data == DATA
assert encoded.parity == PARITY_CONVENTIONAL
encoded = encode_block(DATA, dual_basis=True)
assert encoded.data == DATA
assert encoded.parity == PARITY_DUAL_BASIS
data = [x ^ e for x, e in zip(DATA, DATA_ERROR_MASK)]
parity_conventional = [x ^ e for x, e in zip(PARITY_CONVENTIONAL, PARITY_ERROR_MASK)]
parity_dual_basis = [x ^ e for x, e in zip(PARITY_DUAL_BASIS, PARITY_ERROR_MASK)]
decoded = decode_block(Encoded(data=data, parity=parity_conventional))
assert decoded == DATA
#assert corrected_byte_count == 16 # TODO API
decoded = decode_block(Encoded(data=data, parity=parity_dual_basis), dual_basis=True)
assert decoded == DATA
#assert corrected_byte_count == 16 # TODO API
def _self_bench():
DATA = list(range(223))
DATA_ERROR_MASK = [
0x58,0x00,0xA3,0x00,0x00,0xCD,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x0D,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0xCA,0x96,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x1B,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0xA2,0x00,0xAC,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0xB9,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0xE5,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x94,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0xC3,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x97,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00
]
PARITY_ERROR_MASK = [
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x7A,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x29,0x00,0x00,0x00,0x00,
]
import timeit
def run(stmt, **globals):
timer = timeit.Timer(stmt, globals=globals)
i = 1
time_taken = 0
while time_taken < 2:
for j in 1, 2, 5:
number = i * j
time_taken = timer.timeit(number)
if time_taken >= 2: break
i *= 10
print(stmt, f'[{time_taken / number * 1000:.6f}ms per iter]')
run('encode_block(data)', data=list(range(223)), encode_block=encode_block)
encoded = encode_block(list(range(223)))
run('decode_block(encoded)', encoded=encoded, decode_block=decode_block)
encoded.data = [x ^ e for x, e in zip(encoded.data, DATA_ERROR_MASK)]
encoded.parity = [x ^ e for x, e in zip(encoded.parity, PARITY_ERROR_MASK)]
run('decode_block(corrupt_encoded)', corrupt_encoded=encoded, decode_block=decode_block)
if __name__ == '__main__':
_self_test()
_self_bench()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment