MTH 810 Projects

Luke Baumann April 21st, 2021

*Trying out Jupyter slides. No spell check. Apologies in advance.

Outline

1. Polynomial Class
2. Sieve of Eratosthenes for Irreducible Polynomials
3. Fast Number Theoretic Transform and Inverse
4. QR Code Overview

Polynomial Class

Polynomials in $F_{p^\infty}$

Base Implemented Operations

1. __add__: Addition
2. __neg__: Additive inverse
3. __lshift__: Multiplication by x**n
4. __mul__: Multiplication by integers and polynomials
5. __divmod__: Long Division by integers and polynomials

Derived Operators

1. __sub__: Addition + Additive inverse -> Subtraction
2. __floordiv__: Quotient -> Long Division
3. __mod__: Remainder -> Long Division
4. __pow__: Multiplication -> Exponentiation

Polynomial Class

Because $p$ is prime, for $a \in F_p$, $a^{-i} = a^{p - i - 1}$.

Likewise, if $P(x)$ is an irreducible polynomial, for $a \in F_{p^n} / P(x)$, $a^{-i} = a^{p^n - i - 1}$.

Thus, __pow__ gives the multiplicative inverse.

Utility

1. generate_polynomials: Generates all polynomials up to a give degree
2. generate_monomials: Generate all monomials up to a give degree
3. primitive_element: Finding primitive element given a characteristic polynomial

Sieve of Eratosthenes for Irreducible Polynomials

Rudimentary method to find irreducible polynomials below a certain degree.

Review: Sieve of Eratosthenes for Prime Numbers

Rudimentary method to find prime numbers below a certain number.

def sieve_of_eratosthenes_primes(N):
    # All positive numbers except 1 which is defined to be not irreducible 
    potential_primes = set(range(2, N + 1))

    for number in range(N):
        if number in potential_primes:
            # Remove all multiples of number
            potential_primes -= set(number * i
                                    for i in range(2, ceil(N / number) + 1))

    return potential_primes

N = 100
print('Primes under {}'.format(N), sieve_of_eratosthenes_primes(100))

Primes under 100 {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97}

Sieve of Eratosthenes for Irreducible Polynomials

Number of Irreducible Polynomials

Can confirm the number of irreducible polynomials by Gauss's formula where $p$ is the dictionary size, $n$ is the degree of the monomial, $D$ is the set of all positive divisors of $n$, and $\mu$ is Möbius number (1 if number of factors if even, otherwise, -1).

$$ I_{p^n} = \sum \limits_{d \in D} \mu \left( n / d \right) p^d $$

def I_pn(p, n):
    l = [mobius(d)*p**(n // d) for d in divisors(n)]
    return int(sum(l) // n)

def sieve_of_eratosthenes_irreducible_polynomials(p, N):
    # All monomials except 1*x**0 which is defined to be not irreducible 
    potential_irreducible_monomials  = set(PolynomialNP.generate_monomials(p, 1, N))

    for mono in PolynomialNP.generate_monomials(p, N):
        if mono in potential_irreducible_monomials:
            # Remove all multiples of poly
            potential_irreducible_monomials -= set(
                mono * mono_other
                for mono_other in PolynomialNP.generate_monomials(p, 1, N - mono.degree))

    # Place monomials into a dictionary for easy inspection
    ret = {degree: set() for degree in range(1, N + 1)}
    for mono in potential_irreducible_monomials:
        ret[mono.degree].add(mono)
        
    return ret

p = 3
N = 5
irreducible_polynomials = sieve_of_eratosthenes_irreducible_polynomials(p, N)
for n, ip in irreducible_polynomials.items():
    print('Degree', n)
    print('    I_pn:', I_pn(p, n))
    print('    # Found:', len(ip))
    print()
print('Irreducible Polynomials of Degree 2')
print(irreducible_polynomials[2])

Degree 1
    I_pn: 3
    # Found: 3

Degree 2
    I_pn: 3
    # Found: 3

Degree 3
    I_pn: 8
    # Found: 8

Degree 4
    I_pn: 18
    # Found: 18

Degree 5
    I_pn: 48
    # Found: 48

Irreducible Polynomials of Degree 2
{1*x**2 + 1*x**0, 1*x**2 + 1*x**1 + 2*x**0, 1*x**2 + 2*x**1 + 2*x**0}

Fast Number Theoretic Transform and Inverse

Fast computation of a transfrom applied to polynomials over a finite field.

Why Transform

Certain operators and analysis are easier (simpler, faster) to do in a transformed domain than the original domain.

Review: Discrete Fourier Transforms (DFT)

Invertibly transform a sequence $x_n \in \mathbb{C}$ into another sequence $X_k \in \mathbb{C}$ to perform certain operations easier.

Incomplete List of Benefits

1. Finite differences become multiplication by $1 - e^{-j 2 \pi k / N}$.
2. Cummulative sums become division by $1 - e^{-j 2 \pi k / N}$.
3. Convolution becomes point-wise multiplication.
4. Harmonic analysis is more natural

Review: DFT and IDFT

$$ \begin{align} X_{k} =& \sum \limits _{n=0}^{N-1}x_{n} \left( e^{-j 2\pi /N} \right)^{kn} & X = W^+ x\\ x_{n} =& \frac{1}{K} \sum \limits _{k=0}^{K-1}X_{k} \left( e^{j 2\pi /K} \right)^{nk} & x = \frac{1}{K} W^- X \end{align} $$

Results in a Vandermonde matrix with $\gamma^+ = e^{-j 2\pi / N}$ and $\gamma^- e^{j 2\pi / K}$.

$$W^+=\begin{bmatrix} 1&1&1&1&\cdots &1\\ 1&{\gamma^+} &{\gamma^+} ^{2}&{\gamma^+} ^{3}&\cdots &{\gamma^+} ^{N-1}\\1&{\gamma^+}^{2}&{\gamma^+}^{4}&{\gamma^+}^{6}&\cdots &{\gamma^+} ^{(N-1)2}\\1&{\gamma^+} ^{3}&{\gamma^+} ^{6}&{\gamma^+} ^{9}&\cdots &{\gamma^+}^{(N-1)3}\\\vdots &\vdots &\vdots &\vdots &\ddots &\vdots \\1&{\gamma^+} ^{K-1}&{\gamma^+} ^{2(K-1)}&{\gamma^+} ^{3(K-1)}&\cdots &{\gamma^+}^{(N-1)(K-1)}\end{bmatrix}$$

Similar for $W^-$.

Review: Cooley-Tukey FFT Algorithm

Reexpress DFT of size $N = N_1 N_2$ into $N_1$ DFTs of size $N_2$ recursively, i.e., each DFT of size $N_2 = N_3 N_4$ into $N_3$ DFTs of size $N_4$, etc.

E.g., radix-2 special case where $N_1 = 2$ and $N_2 = N / 2$

$$ X_{k}=\underbrace {\sum \limits _{i=0}^{N_2-1}x_{2n} \left( e^{-j 2\pi / N_2} \right)^{kn}} _{\mathrm {DFT\;of\;even-indexed\;part\;of\;} x_{n}}{} + \left(e^{-j 2 \pi / N} \right)^k \underbrace {\sum \limits _{n=0}^{N_2-1}x_{2n+1} \left( e^{-j 2\pi / N_2} \right)^{kn}} _{\mathrm {DFT\;of\;odd-indexed\;part\;of\;} x_{n}} $$$$ \begin{align} X_{k}=&E_{k}+\left( e^{-j 2\pi / N} \right)^k O_{k}\\ X_{k+N_2}=&E_{k}- \left(e^{-j 2\pi / N} \right)^k O_{k} \end{align} $$

Similar for IFFT.

def fft_radix2(x_n):
    """
    Cooley-Tukey Algorithm for Radix-2.
    
    Performs an FFT recusively until the size is no longer a multiple of 2 and then reverts to a DFT.
    """
    N = x_n.size
    if N < MINIMUM_N:
        return dft(x_n)

    # Check to see if N is divisible by 2. If not, perform a DFT
    try:
        N_1, N_2 = smallest_prime_divisor(N, 2)
    except ValueError:
        return dft(x_n)
    
    gamma = np.exp(-1j * 2 * np.pi / N)
    twiddle = gamma**np.arange(N_2)
    
    even = fft_radix2(x_n[::N_1])
    odd = fft_radix2(x_n[1::N_1])

    twiddle_odd = twiddle * odd

    ret = np.zeros(N, dtype=np.complex128)
    ret[:N_2] = even + twiddle_odd
    ret[N_2:] = even - twiddle_odd

    return ret

def ifft_radix2(X_k):
    """
    Cooley-Tukey Algorithm for Radix-2.
    
    Performs an IFFT recusively until the size is no longer a multiple of 2 and then reverts to a IDFT.
    """
    N = X_k.size
    if N < MINIMUM_N:
        return idft(X_k)
    
    # Check to see if K is divisible by 2. If not, perform a IDFT
    try:
        N_1, N_2 = smallest_prime_divisor(N, 2)
    except ValueError:
        return idft(X_k)

    gamma = np.exp(1j * 2 * np.pi / N)
    twiddle = gamma**np.arange(N_2)

    even = ifft_radix2(X_k[::N_1])
    odd = ifft_radix2(X_k[1::N_1])

    twiddle_odd = twiddle * odd

    ret = np.zeros(N, dtype=np.complex128)
    ret[:N_2] = even + twiddle_odd
    ret[N_2:] = even - twiddle_odd

    return ret / N_1

N = 2**10
x_n = rng.random(N)
X_k = fft_radix2(x_n)
x_n1 = ifft_radix2(X_k)

print('FFT Error:', linalg.norm(X_k - np.fft.fft(x_n)) / linalg.norm(np.fft.fft(x_n)))
print('IFFT Error:', linalg.norm(x_n1 - np.fft.ifft(X_k)) / linalg.norm(np.fft.ifft(X_k)))
print('IFFT(FFT(x_n)) Error:', linalg.norm(x_n - x_n1) / linalg.norm(x_n))

FFT Error: 4.201343978483622e-15
IFFT Error: 4.174086104620431e-15
IFFT(FFT(x_n)) Error: 7.020672682930469e-15

Review: Cooley-Tukey FFT Algorithm

More general case where $N = N_1 N_2$ where $N$ is the number of samples (in both time and frequency) at each level of the recursion.

$$ \begin{align} X_{k_1 N_2 + k_2}=& \sum \limits_{n_1=0}^{N_1 - 1} \left(e^{-j 2 \pi / N} \right)^{n_1 \left( N_2 k_1 + k_2 \right)} \left( \sum \limits _{n_2=0}^{N_2-1}x_{N_1 n_2 + n_1} \left( e^{-j 2\pi / N_2} \right)^{k_2 n_2} \right) \\ x_{n_1 N_2 + n_2}=& \frac{1}{N_1} \sum \limits_{k_1=0}^{N_1 - 1} \left(e^{j 2 \pi / N} \right)^{k_1 \left( N_2 n_1 + n_2 \right)} \left( \frac{N_1}{N} \sum \limits _{k_2=0}^{N_2-1}X_{N_1 k_2 + k_1} \left( e^{j 2\pi / N_2} \right)^{n_2 k_2} \right) \end{align} $$

def fft(x_n):
    N = x_n.size
    if N < MINIMUM_N:
        return dft(x_n)
    
    N_1, N_2 = smallest_prime_divisor(N)

    gamma = np.exp(-1j * 2 * np.pi / N)
    
    parts = [fft(x_n[n_1::N_1]) for n_1 in range(N_1)]
    
    ret = np.zeros(N, dtype=np.complex128)
    for k_1 in range(N_1):
        ret[k_1 * N_2:(k_1 + 1) * N_2] = sum(
            gamma**(n_1 * (N_2 * k_1 + np.arange(N_2))) * part
            for n_1, part in enumerate(parts))

    return ret

def ifft(X_k):
    N = X_k.size
    if N < MINIMUM_N:
        return idft(X_k)
    
    N_1, N_2 = smallest_prime_divisor(N)

    gamma = np.exp(1j * 2 * np.pi / N)
    
    parts = [ifft(X_k[k_1::N_1]) for k_1 in range(N_1)]

    ret = np.zeros(N, dtype=np.complex128)
    for n_1 in range(N_1):
        ret[n_1 * N_2:(n_1 + 1) * N_2] = sum(
            gamma**(k_1 * (N_2 * n_1 + np.arange(N_2))) * part
            for k_1, part in enumerate(parts))

    return ret / N_1

N = 2 * 3 * 5 * 7 * 11
x_n = rng.random(N)
X_k = fft(x_n)
x_n1 = ifft(X_k)

print('FFT Error:', linalg.norm(X_k - np.fft.fft(x_n)) / linalg.norm(np.fft.fft(x_n)))
print('IFFT Error:', linalg.norm(x_n1 - np.fft.ifft(X_k)) / linalg.norm(np.fft.ifft(X_k)))
print('IFFT(FFT(x_n)) Error:', linalg.norm(x_n - x_n1) / linalg.norm(x_n))

FFT Error: 3.632963287721464e-14
IFFT Error: 2.8881688456579974e-14
IFFT(FFT(x_n)) Error: 5.745776431090915e-14

N 840
Radix-2 FFT
    Frequency Domain Error: 1.2417472992382048e-13
    Time Domain Error: 2.830460714125949e-14

FFT
    Frequency Domain Error: 1.2417472992382048e-13
    Time Domain Error: 2.830460714125949e-14

Numpy
    Frequency Domain Error: 1.2606068718914065e-13
    Time Domain Error: 3.415706460783746e-16

Number Theoretic Transform (NTT)

Invertibly transform a sequence $a_j \in F_{p^n} / P(x)$ (or with slight modification, a ring of integers modulo $p$) into another sequence $A_i \in F_{p^n} / P(x)$ to perform certain operations easier.

Incomplete List of Benefits

1. Convolution becomes point-wise multiplication.
   1. Multiplication of polynomials over $F_{p^n} / P(x)$
   2. Multiplication of integer polynomials
   3. Multiplication of very large integers
   4. Division of polynomials over $F_{p^n} / P(x)$
2. Decoding Reed-Solomon (and Reed-Muller) codes fast because their generator/parity check matrices are Vandermonde matrices of primitive elements.

NTT (and INTT) over Finite Fields

1. Field $F_{p^n} / P(x)$ where $p$ is prime and $n \in \mathbb{Z}^+$ and $P(x)$ is an order $n$ irreducible polynomial.
2. Multiplicative group ${F^*}_{p^n} / P(x)$ which is composed of nonzero elements of $F_{p^n} / P(x)$.
3. $d d^\prime = p^n - 1$.
4. $w$ is the primitive element of $F_{p^n} / P(x)$.
5. $a_j$ is a sequence of $d$ elements in $F_{p^n} / P(x)$.
6. $A_i$ is a sequence of $d$ elements in $F_{p^n} / P(x)$.

We use the Polynomial class to get elements of $F_{p^n} / P(x)$ by using modulo arithmetic on polynomials.

$$\begin{matrix} A_i =& \sum \limits _{j = 0}^{d - 1} a_j \left( w^{d^\prime} \right)^{i j} & A = W^+ a & \gamma^+ = w^{d^\prime}\\ a_j =& -d^\prime \sum \limits _{i=0}^{d - 1} A_i \left( w^{-d^\prime} \right)^{j i} & a = -d^\prime W^- A & \gamma^- = w^{-d^\prime} \end{matrix}$$

DFT/NTT Comparison

Discrete Fourier Transform	Number Theoretic Transform
$x \in \mathbb{C}$	$a \in F_{p^n} / P(x)$
$X_k = \sum_{n = 0}^{N - 1} x_n \left( e^{-j 2 \pi / N} \right)^{k n}$	$A_i = \sum \limits _{j = 0}^{d - 1} a_j \left( w^{d^\prime} \right)^{i j}$
$x_n = \frac{1}{N} \sum_{k = 0}^{N - 1} X_k \left( e^{j 2 \pi / N} \right)^{n k}$	$a_j = -d^\prime \sum \limits _{i = 0}^{d - 1} A_i \left( w^{-d^\prime} \right)^{j i}$
$W^\pm$ is a Vandermonde matrix with $\gamma^+ = e^{-j 2 \pi / N}$ and $\gamma^- = e^{j 2 \pi / N}$	$W^\pm$ is a Vandermonde matrix with $\gamma^+ = w^{d^\prime}$ and $\gamma^- = w^{-d^\prime}$
$e^{-j 2 \pi / N}$ is the $N$th root-of-unity	$w$ is the primitive element of ${F^*}_{p^n} / P(x)$
Can be computed fast	Can be computed fast

Cooley-Tukey FNTT Algorithm

Skipping the radix-2 specification, the general case is shown with $d = d_1 d_2$ where $d$ is the number of samples (in both domains) at each level of the recursion.

$$ \begin{align} A_{i_1 d_2 + i_2}=& \sum \limits_{j_1=0}^{d_1 - 1} \left(w^{d^\prime} \right)^{j_1 \left( d_2 i_1 + i_2 \right)} \left( \sum \limits _{j_2=0}^{d_2-1} a_{d_1 j_2 + j_1} \left( w^{d^\prime} \right)^{i_2 j_2} \right) \\ a_{j_1 d_2 + j_2}=& d_1^{-1} \sum \limits_{i_1=0}^{d_1 - 1} \left(w^{-d^\prime} \right)^{i_1 \left( d_2 j_1 + j_2 \right)} \left( -d^\prime d_1 \sum \limits _{i_2=0}^{d_2-1} A_{d_1 i_2 + i_1} \left( w^{-d^\prime} \right)^{j_2 i_2} \right) \end{align} $$

THIS IS THE MOST IMPORTANT RESULT FOR THIS WHOLE PROJECT

def fntt(a_j, p, n, characteristic_polynomial):
    d = len(a_j)
    # Base case, for small d, do naive NTT
    if d < MINIMUM_D:
        return ntt(a_j, p, n, characteristic_polynomial)

    d_prime, remainder = divmod(p**n - 1, d)
    if remainder != 0:
        raise ValueError('d d_prime must equal p**n - 1')

    d_1, d_2 = smallest_prime_divisor(d)

    w = PolynomialNP.primitive_element(p, n, characteristic_poly)
    r = w**d_prime

    parts = [fntt(a_j[j_1::d_1], p, n, characteristic_polynomial)
             for j_1 in range(d_1)]

    A_i = [PolynomialNP(p, n) for i in range(d)]
    for i_1 in range(d_1):
        for i_2 in range(d_2):
            A_i[i_1 * d_2 + i_2] += sum(
                (r**(j_1 * (d_2 * i_1 + i_2)) * part[i_2]
                for j_1, part in enumerate(parts)), PolynomialNP(p, n))

    for i in range(d):
        A_i[i] %= characteristic_poly

    return A_i

def fintt(A_i, p, n, characteristic_polynomial):
    d = len(A_i)
    # Base case, for small d, do naive INTT
    if d < MINIMUM_D:
        return intt(A_i, p, n, characteristic_polynomial)

    d_prime, remainder = divmod(p**n - 1, d)
    if remainder != 0:
        raise ValueError('d d_prime must equal p**n - 1')

    d_1, d_2 = smallest_prime_divisor(d)

    w = PolynomialNP.primitive_element(p, n, characteristic_poly)
    r = w**(p**n - d_prime - 1)

    parts = [fintt(A_i[i_1::d_1], p, n, characteristic_polynomial)
             for i_1 in range(d_1)]

    a_j = [PolynomialNP(p, n) for j in range(d)]
    for j_1 in range(d_1):
        for j_2 in range(d_2):
            a_j[j_1 * d_2 + j_2] += sum(
                (r**(i_1 * (d_2 * j_1 + j_2)) * part[j_2]
                for i_1, part in enumerate(parts)), PolynomialNP(p, n))

    for j in range(d):
        a_j[j] = d_1**(p - 1 - 1) * a_j[j] % characteristic_poly

    return a_j

p = 3
n = 2
d_prime = 2
d, remainder = divmod(p**n - 1, d_prime)
if remainder != 0:
    raise ValueError('d d_prime must equal p**n - 1')
characteristic_poly = PolynomialNP(p, [2, 1, 1])
primitive_element = PolynomialNP.primitive_element(p, n, characteristic_poly)

a_j = [PolynomialNP(p, rng.integers(0, p, n)) for i in range(d)]
A_i = fntt(a_j, p, n, characteristic_poly)
a_j1 = fintt(A_i, p, n, characteristic_poly)

print('FNTT Error:', [Ai - Ai1 for Ai, Ai1 in zip(A_i, ntt(a_j, p, n, characteristic_poly))])
print('FINTT Error:', [aj - aj1 for aj, aj1 in zip(a_j, intt(A_i, p, n, characteristic_poly))])
print('FINTT(FNTT) Error:', [aj - aj1 for aj, aj1 in zip(a_j, a_j1)])

FNTT Error: [0, 0, 0, 0]
FINTT Error: [0, 0, 0, 0]
FINTT(FNTT) Error: [0, 0, 0, 0]

My Open Questions

1. What is the analog to finite difference? $A_j \left(1 - \left( w^{d^\prime} \right)^j \right)$
2. What is the analog to cummulative sum? $A_j \left(1 - \left( w^{d^\prime} \right)^j \right)^{-1}$
3. What is the analog to harmonic analysis? The interpretation of $A_j$.

QR Code Overview

Most of this is taken from Wikipedia but a very thorough tutorial is given here: https://www.thonky.com/qr-code-tutorial/

QR Code key facts

Quick Response Code is 2D barcode (machine-readable, optical label) invented in 1994.

• Micro QR Code sizes: 11x11 to 17x17
• QR Code sizes: 21x21 to 177x177
• Can encode numeric data (0-9), alphanumeric (0-9, a-z, A-Z, " $%*+-./:", byte data, and Kanji characters)
• Four levels of Reed-Solomon error correction (L: 7%, M: 15%, Q: 25%, H: 30%)
• Max Micro QR Code data: 35 numeric, 21 alphanumeric, 15 bytes, 9 Kanji.
• Max QR Code data: 7089 numeric, 4296 alphanumeric, 2953 bytes, 1817 Kanji.
• Reflection, rotation, color inversion independent
• Introduce intentional errors to incorporate colors, logos, art, etc.

FYI: QR Codes flip the terminology for blocks and codewords.

The number of erasures and errors that can be corrected is given by $e + 2 t \le d - p$ where $e$ is the number of erasures, $t$ is the number of errors, $d$ is the number of error correcting codewords, and $p$ is the number of misdecode protection codewords. Not covered in our course, $p$ is the number of extra codewords that can help detect errors but not correct them, $n - k = d + p$.

Version	Total number of codewords	Error correction Level	Number of error correction codewords	Value of $p$	Number of error correction blocks	Error correction code</br>per block $(n, k, d)_{2^8}$
M1	5	Error detection only	2	2	1	(5, 3, 0)
M4	24	L M Q	8 10 14	2 0 0	1	(24, 16, 3)a (24, 14, 5) (24, 10, 7)
1	26	L M Q H	7 10 13 17	3 2 1 1	1	(26, 19, 2)a (26, 16, 4)a (26, 13, 6)a (26, 9, 8)a
10	346	L M Q H	72 130 192 224	0	2 2 4 1 6 2 6 2	(86, 68, 9) (87, 69, 9) (69, 43, 13) (70, 44, 13) (43, 19, 12) (44, 20, 12) (43, 15, 14) (44, 16, 14)
20	1085	L M Q H	224 416 600 700	0	3 5 3 13 15 5 15 10	(135, 107, 14) (136, 108, 14) (67, 41, 13) (68, 42, 13) (54, 24, 15) (55, 25, 15) (43, 15, 14) (44, 16, 14)
40	3706	L M Q H	750 1372 2040 2430	0	19 6 18 31 34 34 20 61	(148, 118, 15) (149, 119, 15) (75, 47, 14) (76, 48, 14) (54, 24, 15) (55, 25, 15) (45, 15, 15) (46, 16, 15)

a: Error correction capacity is less than half the number of error correction codewords to reduce the probability of misdecodes.

Thanks!