r"""
Peterson and Brown's cyclic redundancy check (1961)
===================================================

W. Wesley Peterson and D. T. Brown showed how to detect errors with a few
shift-register stages and exclusive-OR gates. Treat the message as a
polynomial over GF(2), multiply by :math:`x^r`, divide by a generator
:math:`g(x)` of degree :math:`r`, and send the remainder after the message
(:func:`~engineeringkit.communications.crc_encode`). The receiver divides
again and declares an error if anything is left
(:func:`~engineeringkit.communications.crc_check`). Every burst of :math:`r`
or fewer errors is caught, and a random error pattern slips through with
probability :math:`2^{-r}`. CRCs now guard Ethernet frames, disks and USB
packets.
"""

# %%
import matplotlib.pyplot as plt
import numpy as np

from engineeringkit.communications import CRC16_CCITT, crc_check, crc_encode, crc_remainder

# %%
# Long division by x^3 + x + 1
# ----------------------------

message = [1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0]
print("remainder:", crc_remainder(message, 0b1011).tolist())
assert crc_remainder(message, 0b1011).tolist() == [1, 0, 0]

# The catalogue check value of CRC-16/XMODEM: the ASCII string "123456789".
check = crc_remainder(np.unpackbits(np.frombuffer(b"123456789", dtype=np.uint8)), CRC16_CCITT)
print(f"CRC-16-CCITT of '123456789' = 0x{int(''.join(map(str, check)), 2):04X}")
assert int("".join(map(str, check)), 2) == 0x31C3

# %%
# Bursts are always caught; random patterns almost always
# -------------------------------------------------------

rng = np.random.default_rng(1961)
frame = crc_encode(rng.integers(0, 2, 96), CRC16_CCITT)
missed_bursts = 0
for length in range(1, 17):
    for start in range(frame.size - length + 1):
        pattern = np.zeros(frame.size, dtype=int)
        pattern[start] = pattern[start + length - 1] = 1
        pattern[start + 1 : start + length - 1] = rng.integers(0, 2, max(length - 2, 0))
        missed_bursts += crc_check(frame ^ pattern, CRC16_CCITT)
print("undetected bursts of length <= 16:", missed_bursts)
assert missed_bursts == 0

short = crc_encode(rng.integers(0, 2, 12), 0b10011)  # CRC-4, so misses are frequent enough to count
trials = 40_000
missed = sum(crc_check(short ^ rng.integers(0, 2, short.size), 0b10011) for _ in range(trials))
print(f"random errors missed by CRC-4: {missed / trials:.4f} (2^-4 = {2**-4:.4f})")
assert abs(missed / trials - 2**-4) < 0.01

fig, ax = plt.subplots(figsize=(6, 3.6))
r = np.arange(1, 33)
ax.semilogy(r, 2.0**-r, label="random pattern undetected, $2^{-r}$")
ax.plot([4], [missed / trials], "o", label="CRC-4, simulated")
ax.set(xlabel="check bits r", ylabel="probability")
ax.legend()
fig.tight_layout()

plt.show()