Shared2018-07-12-215555.sagewsOpen in CoCalc
McEliece Toy Implementation
# ----------------------------------
# TOY McEliece (code-based) ENCRYPTION SCHEME
# This is purely to demonstrate my understanding of code-based crypto
#
# I used the specification for keygen, enc and dec provided in the following wikipedia page and in the original McEliece specification paper:
# https://en.wikipedia.org/wiki/McEliece_cryptosystem
# https://ipnpr.jpl.nasa.gov/progress_report2/42-44/44N.PDF
#
# I also used the version of Patterson's algorithm for binary goppa decoding here (PAGE 5):
# https://pdfs.semanticscholar.org/2754/dd1eef3f1b50b88df0b6689dc0478d22c34e.pdf
# ----------------------------------

# Sage functions for converting from ascii to binary strings and vice versa
from sage.crypto.util import bin_to_ascii, ascii_to_bin
import random

#################################################
# 1 - Parameters
# Parameters for system
n = 31; k = 16; t = 3
#n = 1024; k = 524; t = 50 # McEliece's original parameter set
#n = 6960; k = 5413; t = 119 # Classic McEliece parameters for PQ security of 2^128 bits
m = (n - k) / t

F.<x> = GF(2) # Used for computing binary matrices e.g. S and P
phi.<x> = GF(2^m) # Used for computing the binary goppa code
sequence = [x^i for i in range(n)] # part of the goppa code along with the polynomial g
PolyRing = PolynomialRing(phi, 'z')
z = PolyRing.gen()

#################################################
# 2 - Error polynomial of weight t - to be added to the encoded message as part of encryption
def error_vector():
    error_vector = vector(F, n)
    #return error_vector
    points = random.sample(range(n), t)

    for i in range(t):
        error_vector[points[i]] += 1

    return error_vector

#################################################
# 3 - Encoding a message into a polynomial
def msg_encode(F, msg):
    S = ascii_to_bin(msg)
    assert len(str(S)) < k
    msg_lngth = len(str(S))

    msg = vector(F, k)

    for i in range(len(str(S))):
        msg[i] = int(str(S[i]))

    return msg, msg_lngth

#################################################
# 4 - Decoding algorithm which dec() uses to obtain the original message
def msg_decode(msg, l):
    return bin_to_ascii(''.join(str(msg[0, i]) for i in range(l)))

#################################################
# 5 - Keygen (outputs private key and public key)

# Generate binary goppa polynomial g(x)
def goppa_poly():

    # Randomly chosen Monic irreducible polynomial of degree t
    return phi[z].irreducible_element(t, algorithm = "random")

# Generate k x k binary invertible matrix
def non_singular_matrix(F, k):

    # The code generates a random matrix and if it is not invertible (i.e. is singular) then a new matrix is generated until it is
    while 1:
        try:
            r = random_matrix(F, k, k)
            if r.inverse() != ZeroDivisionError:
                break
        except ZeroDivisionError:
            continue
    return r

# Generate n x n random permutation matrix (without this the scheme would reveal parts of the plaintext)
def rand_perm(F, n):
    r = matrix(F, n, n)

    # Sample a unique set of randomly chosen elements and use these to set the corresponding row element to 1
    points = random.sample(range(n), n)

    for i in range(n):
        r[i, points[i]] = 1

    # Make sure r is a random permutation matrix before returning it
    v = vector(F, n)
    for i in range(n):
        v += r[i]
    # Ensure r is correctly computed
    assert v + v == 0

    return r

def keygen(g):

    # Build the Vandermonde matrix - part of the parity check matrix
    H_g = matrix([[sequence[j] ^ (i) for j in range(n)] for i in range(t)]);

    # Multiply by the diagonal matrix of inverted elements of the goppa code
    H_g *= diagonal_matrix([1/g(sequence[i]) for i in range(n)])

    H_Goppa = matrix(F, m * H_g.nrows(), H_g.ncols())

    for i in range(H_g.nrows()):
        for j in range(H_g.ncols()):

            # For each matrix element, obtain the bitstring representation of the polynomial
            # Then reverse the order and add a set number of 0's to the end
            # Finally convert the result to a list and then a vector
            bitstring = bin(H_g[i, j].integer_representation())[2:];
            bitstring = bitstring[::-1];
            bitstring = bitstring + "".join(str(0) for i in range(m - len(bitstring)))
            bitstring = list(bitstring);

            H_Goppa[m * i:m * (i + 1), j] = vector(map(int, bitstring))

    P = rand_perm(F, n)
    S = non_singular_matrix(F, k)
    G = H_Goppa.right_kernel().basis_matrix()
    G_hat = S * G * P

    return (G_hat, t), (S, G, P), H_Goppa

#################################################
# 6 - Encryption function - multiply the message vector by the public key and add the error vector
def enc(pk, msg, e):
    return (msg * pk[0]) + e

#################################################
# 7 - Decryption
# There are three stages to this function:
# - Apply the inverse of the permutation matrix P
# - Apply Patterson's decoding algorithm to recover and return msg'

def split(g):
    Phi = g.parent()
    g0 = Phi([sqrt(i) for i in g.list()[0::2]]) # Choose terms in even positions in g
    g1 = Phi([sqrt(i) for i in g.list()[1::2]]) # Choose terms in odd positions in g

    # Ensure the splitting worked
    assert g0^2 + z*(g1^2) == g

    return (g0, g1)

def dec(sk, ct, g, H):

    # Remove the random permutation P
    ct *= sk[2].inverse()

    # METHOD 1 SYNDROME COMPUTATION
    syndrome = H * ct.column()
    syn_poly = PolyRing(0)

    for i in range(syndrome.ncols()):
        syn_poly += syndrome[i, 0]*(z^i)

    #syn_poly = syn_poly.mod(g)

    # METHOD 2 SYNDROME COMPUTATION - Both methods produce the same result
    #syn = vector(g.parent(), len(sequence))
    #for i in range(len(sequence)):
    #    syn[i] = ((z - sequence[i]).inverse_mod(g))
    #
    #syn *= ct.column()
    #syn_poly = PolyRing(0)
    #for i in range(len(syn)):
    #    syn_poly += syn[i]*(z^i)
    #
    #syn_poly = syn_poly.mod(g)

    if syn_poly == 0:
        print "no errors added"
        return ct

    T = syn_poly.inverse_mod(g) # When no errors are added to the codeword this if statement should be triggered, but it isnt so something is either wrong with the syndrome or the goppa polynomial
    (T0, T1) = split(T + z)
    assert T0^2 + z*(T1^2) == T + z
    (g0, g1) = split(g)
    w = g0 * (g1.inverse_mod(g))
    R = (T0 + (w * T1)).mod(g)
    (d, u, v) = xgcd(PolyRing(1), PolyRing(R.list()))
    a = g * u
    b = g * v
    if a == (b * R).mod(g):
        sigma = (a ^ 2) + (z * (b ^ 2))
    else:
        print "error"
        sigma = (a ^ 2) + (z * (b ^ 2))

    for i in range(len(sequence)):
        if (sigma(sequence[i]) == 0):
            print "Error in position: " + str(i)
            ct[i] += 1

    return ct

#################################################
# 8 - Main

# Compute the Goppa code (as a polynomial)
g = goppa_poly()

# Compute the secret key and public key tuples
pk, sk, H = keygen(g)

# Compute the error vector of weight t
e = error_vector()

# Encode a message into a vector of length k
message = 'a'; print "String to be encrypted: " + message
#message = 'McEliece does not work!'; print "String to be encrypted: " + message
msg, msg_lngth = msg_encode(F, message); print msg

# Encrypt the message using the public key
ct = enc(pk, msg, e)

# Decode the ciphertext message using the Patterson algorithm
msg2 = dec(sk, ct, g, H)
msg2 = (sk[1].solve_left(msg2))*sk[0].inverse(); print msg2

# Decode the message vector back into a string
dec_msg = msg_decode(msg2, msg_lngth); print "Decrypted String: " + dec_msg
String to be encrypted: a (0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0) error
Error in lines 106-106 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1055, in execute exec compile(block+'\n', '', 'single', flags=compile_flags) in namespace, locals File "", line 1, in <module> File "sage/matrix/matrix2.pyx", line 205, in sage.matrix.matrix2.Matrix.solve_left (build/cythonized/sage/matrix/matrix2.c:5822) raise ValueError(str(e).replace('row', 'column')) ValueError: matrix equation has no solutions