 diff --git a/CMakeLists.txt b/CMakeLists.txt
-index 7852b3f..6120cb8 100644
+index 7852b3f..03ba98f 100644
 --- a/CMakeLists.txt
 +++ b/CMakeLists.txt
-@@ -2,6 +2,7 @@ cmake_minimum_required(VERSION 3.5)
+@@ -2,12 +2,15 @@ cmake_minimum_required(VERSION 3.5)
  project(wirehair)
  
  set(CMAKE_CXX_STANDARD 11)

 
 add_executable(unit_test ${UNIT_TEST_SOURCE_FILES})
 target_link_libraries(unit_test wirehair)
@@ -88,9 +110,14 @@ add_executable(gen_tables ${GEN_TABLES})
 
 include(GNUInstallDirs)
 

+install(EXPORT ${targetsFile} DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake ${nameSpace})
+install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+install(DIRECTORY python DESTINATION ${CMAKE_INSTALL_PREFIX} ${optionalInstall})
diff --git a/WirehairTools.h b/WirehairTools.h
index 8340df5..f8b3251 100644
--- a/WirehairTools.h
+++ b/WirehairTools.h
@@ -31,7 +31,7 @@
 #define WIREHAIR_TOOLS_H
 
 #include <wirehair/wirehair.h>
-#include "gf256.h"
+#include <wirehair/gf256.h>
 #include <new> // std::nothrow
 
 // Compiler-specific debug break
diff --git a/cm256.cpp b/cm256.cpp
new file mode 100644
index 0000000..b51ff8f
--- /dev/null
+++ b/cm256.cpp
@@ -0,0 +1,603 @@
+/*
+	Copyright (c) 2015 Christopher A. Taylor.  All rights reserved.
+
+	Redistribution and use in source and binary forms, with or without
+	modification, are permitted provided that the following conditions are met:
+
+	* Redistributions of source code must retain the above copyright notice,
+	  this list of conditions and the following disclaimer.
+	* Redistributions in binary form must reproduce the above copyright notice,
+	  this list of conditions and the following disclaimer in the documentation
+	  and/or other materials provided with the distribution.
+	* Neither the name of CM256 nor the names of its contributors may be
+	  used to endorse or promote products derived from this software without
+	  specific prior written permission.
+
+	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+	AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+	ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+	LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+	SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+	INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+	CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+	ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+	POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#include <wirehair/cm256.h>
+
+
+/*
+    GF(256) Cauchy Matrix Overview
+
+    As described on Wikipedia, each element of a normal Cauchy matrix is defined as:
+
+        a_ij = 1 / (x_i - y_j)
+        The arrays x_i and y_j are vector parameters of the matrix.
+        The values in x_i cannot be reused in y_j.
+
+    Moving beyond the Wikipedia...
+
+    (1) Number of rows (R) is the range of i, and number of columns (C) is the range of j.
+
+    (2) Being able to select x_i and y_j makes Cauchy matrices more flexible in practice
+        than Vandermonde matrices, which only have one parameter per row.
+
+    (3) Cauchy matrices are always invertible, AKA always full rank, AKA when treated as
+        as linear system y = M*x, the linear system has a single solution.
+
+    (4) A Cauchy matrix concatenated below a square CxC identity matrix always has rank C,
+        Meaning that any R rows can be eliminated from the concatenated matrix and the
+        matrix will still be invertible.  This is how Reed-Solomon erasure codes work.
+
+    (5) Any row or column can be multiplied by non-zero values, and the resulting matrix
+        is still full rank.  This is true for any matrix, since it is effectively the same
+        as pre and post multiplying by diagonal matrices, which are always invertible.
+
+    (6) Matrix elements with a value of 1 are much faster to operate on than other values.
+        For instance a matrix of [1, 1, 1, 1, 1] is invertible and much faster for various
+        purposes than [2, 2, 2, 2, 2].
+
+    (7) For GF(256) matrices, the symbols in x_i and y_j are selected from the numbers
+        0...255, and so the number of rows + number of columns may not exceed 256.
+        Note that values in x_i and y_j may not be reused as stated above.
+
+    In summary, Cauchy matrices
+        are preferred over Vandermonde matrices.  (2)
+        are great for MDS erasure codes.  (3) and (4)
+        should be optimized to include more 1 elements.  (5) and (6)
+        have a limited size in GF(256), rows+cols <= 256.  (7)
+*/
+
+
+//-----------------------------------------------------------------------------
+// Initialization
+
+extern "C" int cm256_init_(int version)
+{
+    if (version != CM256_VERSION)
+    {
+        // User's header does not match library version
+        return -10;
+    }
+
+    // Return error code from GF(256) init if required
+    return gf256_init();
+}
+
+
+/*
+    Selected Cauchy Matrix Form
+
+    The matrix consists of elements a_ij, where i = row, j = column.
+    a_ij = 1 / (x_i - y_j), where x_i and y_j are sets of GF(256) values
+    that do not intersect.
+
+    We select x_i and y_j to just be incrementing numbers for the
+    purposes of this library.  Further optimizations may yield matrices
+    with more 1 elements, but the benefit seems relatively small.
+
+    The x_i values range from 0...(originalCount - 1).
+    The y_j values range from originalCount...(originalCount + recoveryCount - 1).
+
+    We then improve the Cauchy matrix by dividing each column by the
+    first row element of that column.  The result is an invertible
+    matrix that has all 1 elements in the first row.  This is equivalent
+    to a rotated Vandermonde matrix, so we could have used one of those.
+
+    The advantage of doing this is that operations involving the first
+    row will be extremely fast (just memory XOR), so the decoder can
+    be optimized to take advantage of the shortcut when the first
+    recovery row can be used.
+
+    First row element of Cauchy matrix for each column:
+    a_0j = 1 / (x_0 - y_j) = 1 / (x_0 - y_j)
+
+    Our Cauchy matrix sets first row to ones, so:
+    a_ij = (1 / (x_i - y_j)) / a_0j
+    a_ij = (y_j - x_0) / (x_i - y_j)
+    a_ij = (y_j + x_0) div (x_i + y_j) in GF(256)
+*/
+
+// This function generates each matrix element based on x_i, x_0, y_j
+// Note that for x_i == x_0, this will return 1, so it is better to unroll out the first row.
+static GF256_FORCE_INLINE unsigned char GetMatrixElement(unsigned char x_i, unsigned char x_0, unsigned char y_j)
+{
+    return gf256_div(gf256_add(y_j, x_0), gf256_add(x_i, y_j));
+}
+
+
+//-----------------------------------------------------------------------------
+// Encoding
+
+extern "C" void cm256_encode_block(
+    cm256_encoder_params params, // Encoder parameters
+    cm256_block* originals,      // Array of pointers to original blocks
+    int recoveryBlockIndex,      // Return value from cm256_get_recovery_block_index()
+    void* recoveryBlock)         // Output recovery block
+{
+    // If only one block of input data,
+    if (params.OriginalCount == 1)
+    {
+        // No meaningful operation here, degenerate to outputting the same data each time.
+
+        memcpy(recoveryBlock, originals[0].Block, params.BlockBytes);
+        return;
+    }
+    // else OriginalCount >= 2:
+
+    // Unroll first row of recovery matrix:
+    // The matrix we generate for the first row is all ones,
+    // so it is merely a parity of the original data.
+    if (recoveryBlockIndex == params.OriginalCount)
+    {
+        gf256_addset_mem(recoveryBlock, originals[0].Block, originals[1].Block, params.BlockBytes);
+        for (int j = 2; j < params.OriginalCount; ++j)
+        {
+            gf256_add_mem(recoveryBlock, originals[j].Block, params.BlockBytes);
+        }
+        return;
+    }
+
+    // TBD: Faster algorithms seem to exist for computing this matrix-vector product.
+
+    // Start the x_0 values arbitrarily from the original count.
+    const uint8_t x_0 = static_cast<uint8_t>(params.OriginalCount);
+
+    // For other rows:
+    {
+        const uint8_t x_i = static_cast<uint8_t>(recoveryBlockIndex);
+
+        // Unroll first operation for speed
+        {
+            const uint8_t y_0 = 0;
+            const uint8_t matrixElement = GetMatrixElement(x_i, x_0, y_0);
+
+            gf256_mul_mem(recoveryBlock, originals[0].Block, matrixElement, params.BlockBytes);
+        }
+
+        // For each original data column,
+        for (int j = 1; j < params.OriginalCount; ++j)
+        {
+            const uint8_t y_j = static_cast<uint8_t>(j);
+            const uint8_t matrixElement = GetMatrixElement(x_i, x_0, y_j);
+
+            gf256_muladd_mem(recoveryBlock, matrixElement, originals[j].Block, params.BlockBytes);
+        }
+    }
+}
+
+extern "C" int cm256_encode(
+    cm256_encoder_params params, // Encoder params
+    cm256_block* originals,      // Array of pointers to original blocks
+    void* recoveryBlocks)        // Output recovery blocks end-to-end
+{
+    // Validate input:
+    if (params.OriginalCount <= 0 ||
+        params.RecoveryCount <= 0 ||
+        params.BlockBytes <= 0)
+    {
+        return -1;
+    }
+    if (params.OriginalCount + params.RecoveryCount > 256)
+    {
+        return -2;
+    }
+    if (!originals || !recoveryBlocks)
+    {
+        return -3;
+    }
+
+    uint8_t* recoveryBlock = static_cast<uint8_t*>(recoveryBlocks);
+
+    for (int block = 0; block < params.RecoveryCount; ++block, recoveryBlock += params.BlockBytes)
+    {
+        cm256_encode_block(params, originals, (params.OriginalCount + block), recoveryBlock);
+    }
+
+    return 0;
+}
+
+
+//-----------------------------------------------------------------------------
+// Decoding
+
+struct CM256Decoder
+{
+    // Encode parameters
+    cm256_encoder_params Params;
+
+    // Recovery blocks
+    cm256_block* Recovery[256];
+    int RecoveryCount;
+
+    // Original blocks
+    cm256_block* Original[256];
+    int OriginalCount;
+
+    // Row indices that were erased
+    uint8_t ErasuresIndices[256];
+
+    // Initialize the decoder
+    bool Initialize(cm256_encoder_params& params, cm256_block* blocks);
+
+    // Decode m=1 case
+    void DecodeM1();
+
+    // Decode for m>1 case
+    void Decode();
+
+    // Generate the LU decomposition of the matrix
+    void GenerateLDUDecomposition(uint8_t* matrix_L, uint8_t* diag_D, uint8_t* matrix_U);
+};
+
+bool CM256Decoder::Initialize(cm256_encoder_params& params, cm256_block* blocks)
+{
+    Params = params;
+
+    cm256_block* block = blocks;
+    OriginalCount = 0;
+    RecoveryCount = 0;
+
+    // Initialize erasures to zeros
+    for (int ii = 0; ii < params.OriginalCount; ++ii)
+    {
+        ErasuresIndices[ii] = 0;
+    }
+
+    // For each input block,
+    for (int ii = 0; ii < params.OriginalCount; ++ii, ++block)
+    {
+        int row = block->Index;
+
+        // If it is an original block,
+        if (row < params.OriginalCount)
+        {
+            Original[OriginalCount++] = block;
+
+            if (ErasuresIndices[row] != 0)
+            {
+                // Error out if two row indices repeat
+                return false;
+            }
+
+            ErasuresIndices[row] = 1;
+        }
+        else
+        {
+            Recovery[RecoveryCount++] = block;
+        }
+    }
+
+    // Identify erasures
+    for (int ii = 0, indexCount = 0; ii < 256; ++ii)
+    {
+        if (!ErasuresIndices[ii])
+        {
+            ErasuresIndices[indexCount] = static_cast<uint8_t>( ii );
+
+            if (++indexCount >= RecoveryCount)
+            {
+                break;
+            }
+        }
+    }
+
+    return true;
+}
+
+void CM256Decoder::DecodeM1()
+{
+    // XOR all other blocks into the recovery block
+    uint8_t* outBlock = static_cast<uint8_t*>(Recovery[0]->Block);
+    const uint8_t* inBlock = nullptr;
+
+    // For each block,
+    for (int ii = 0; ii < OriginalCount; ++ii)
+    {
+        const uint8_t* inBlock2 = static_cast<const uint8_t*>(Original[ii]->Block);
+
+        if (!inBlock)
+        {
+            inBlock = inBlock2;
+        }
+        else
+        {
+            // outBlock ^= inBlock ^ inBlock2
+            gf256_add2_mem(outBlock, inBlock, inBlock2, Params.BlockBytes);
+            inBlock = nullptr;
+        }
+    }
+
+    // Complete XORs
+    if (inBlock)
+    {
+        gf256_add_mem(outBlock, inBlock, Params.BlockBytes);
+    }
+
+    // Recover the index it corresponds to
+    Recovery[0]->Index = ErasuresIndices[0];
+}
+
+// Generate the LU decomposition of the matrix
+void CM256Decoder::GenerateLDUDecomposition(uint8_t* matrix_L, uint8_t* diag_D, uint8_t* matrix_U)
+{
+    // Schur-type-direct-Cauchy algorithm 2.5 from
+    // "Pivoting and Backward Stability of Fast Algorithms for Solving Cauchy Linear Equations"
+    // T. Boros, T. Kailath, V. Olshevsky
+    // Modified for practical use.  I folded the diagonal parts of U/L matrices into the
+    // diagonal one to reduce the number of multiplications to perform against the input data,
+    // and organized the triangle matrices in memory to allow for faster SSE3 GF multiplications.
+
+    // Matrix size NxN
+    const int N = RecoveryCount;
+
+    // Generators
+    uint8_t g[256], b[256];
+    for (int i = 0; i < N; ++i)
+    {
+        g[i] = 1;
+        b[i] = 1;
+    }
+
+    // Temporary buffer for rotated row of U matrix
+    // This allows for faster GF bulk multiplication
+    uint8_t rotated_row_U[256];
+    uint8_t* last_U = matrix_U + ((N - 1) * N) / 2 - 1;
+    int firstOffset_U = 0;
+
+    // Start the x_0 values arbitrarily from the original count.
+    const uint8_t x_0 = static_cast<uint8_t>(Params.OriginalCount);
+
+    // Unrolling k = 0 just makes it slower for some reason.
+    for (int k = 0; k < N - 1; ++k)
+    {
+        const uint8_t x_k = Recovery[k]->Index;
+        const uint8_t y_k = ErasuresIndices[k];
+
+        // D_kk = (x_k + y_k)
+        // L_kk = g[k] / (x_k + y_k)
+        // U_kk = b[k] * (x_0 + y_k) / (x_k + y_k)
+        const uint8_t D_kk = gf256_add(x_k, y_k);
+        const uint8_t L_kk = gf256_div(g[k], D_kk);
+        const uint8_t U_kk = gf256_mul(gf256_div(b[k], D_kk), gf256_add(x_0, y_k));
+
+        // diag_D[k] = D_kk * L_kk * U_kk
+        diag_D[k] = gf256_mul(D_kk, gf256_mul(L_kk, U_kk));
+
+        // Computing the k-th row of L and U
+        uint8_t* row_L = matrix_L;
+        uint8_t* row_U = rotated_row_U;
+        for (int j = k + 1; j < N; ++j)
+        {
+            const uint8_t x_j = Recovery[j]->Index;
+            const uint8_t y_j = ErasuresIndices[j];
+
+            // L_jk = g[j] / (x_j + y_k)
+            // U_kj = b[j] / (x_k + y_j)
+            const uint8_t L_jk = gf256_div(g[j], gf256_add(x_j, y_k));
+            const uint8_t U_kj = gf256_div(b[j], gf256_add(x_k, y_j));
+
+            *matrix_L++ = L_jk;
+            *row_U++ = U_kj;
+
+            // g[j] = g[j] * (x_j + x_k) / (x_j + y_k)
+            // b[j] = b[j] * (y_j + y_k) / (y_j + x_k)
+            g[j] = gf256_mul(g[j], gf256_div(gf256_add(x_j, x_k), gf256_add(x_j, y_k)));
+            b[j] = gf256_mul(b[j], gf256_div(gf256_add(y_j, y_k), gf256_add(y_j, x_k)));
+        }
+
+        // Do these row/column divisions in bulk for speed.
+        // L_jk /= L_kk
+        // U_kj /= U_kk
+        const int count = N - (k + 1);
+        gf256_div_mem(row_L, row_L, L_kk, count);
+        gf256_div_mem(rotated_row_U, rotated_row_U, U_kk, count);
+
+        // Copy U matrix row into place in memory.
+        uint8_t* output_U = last_U + firstOffset_U;
+        row_U = rotated_row_U;
+        for (int j = k + 1; j < N; ++j)
+        {
+            *output_U = *row_U++;
+            output_U -= j;
+        }
+        firstOffset_U -= k + 2;
+    }
+
+    // Multiply diagonal matrix into U
+    uint8_t* row_U = matrix_U;
+    for (int j = N - 1; j > 0; --j)
+    {
+        const uint8_t y_j = ErasuresIndices[j];
+        const int count = j;
+
+        gf256_mul_mem(row_U, row_U, gf256_add(x_0, y_j), count);
+        row_U += count;
+    }
+
+    const uint8_t x_n = Recovery[N - 1]->Index;
+    const uint8_t y_n = ErasuresIndices[N - 1];
+
+    // D_nn = 1 / (x_n + y_n)
+    // L_nn = g[N-1]
+    // U_nn = b[N-1] * (x_0 + y_n)
+    const uint8_t L_nn = g[N - 1];
+    const uint8_t U_nn = gf256_mul(b[N - 1], gf256_add(x_0, y_n));
+
+    // diag_D[N-1] = L_nn * D_nn * U_nn
+    diag_D[N - 1] = gf256_div(gf256_mul(L_nn, U_nn), gf256_add(x_n, y_n));
+}
+
+void CM256Decoder::Decode()
+{
+    // Matrix size is NxN, where N is the number of recovery blocks used.
+    const int N = RecoveryCount;
+
+    // Start the x_0 values arbitrarily from the original count.
+    const uint8_t x_0 = static_cast<uint8_t>(Params.OriginalCount);
+
+    // Eliminate original data from the the recovery rows
+    for (int originalIndex = 0; originalIndex < OriginalCount; ++originalIndex)
+    {
+        const uint8_t* inBlock = static_cast<const uint8_t*>(Original[originalIndex]->Block);
+        const uint8_t inRow = Original[originalIndex]->Index;
+
+        for (int recoveryIndex = 0; recoveryIndex < N; ++recoveryIndex)
+        {
+            uint8_t* outBlock = static_cast<uint8_t*>(Recovery[recoveryIndex]->Block);
+            const uint8_t x_i = Recovery[recoveryIndex]->Index;
+            const uint8_t y_j = inRow;
+            const uint8_t matrixElement = GetMatrixElement(x_i, x_0, y_j);
+
+            gf256_muladd_mem(outBlock, matrixElement, inBlock, Params.BlockBytes);
+        }
+    }
+
+    // Allocate matrix
+    static const int StackAllocSize = 2048;
+    uint8_t stackMatrix[StackAllocSize];
+    uint8_t* dynamicMatrix = nullptr;
+    uint8_t* matrix = stackMatrix;
+    const int requiredSpace = N * N;
+    if (requiredSpace > StackAllocSize)
+    {
+        dynamicMatrix = new uint8_t[requiredSpace];
+        matrix = dynamicMatrix;
+    }
+
+    /*
+        Compute matrix decomposition:
+
+            G = L * D * U
+
+        L is lower-triangular, diagonal is all ones.
+        D is a diagonal matrix.
+        U is upper-triangular, diagonal is all ones.
+    */
+    uint8_t* matrix_U = matrix;
+    uint8_t* diag_D = matrix_U + (N - 1) * N / 2;
+    uint8_t* matrix_L = diag_D + N;
+    GenerateLDUDecomposition(matrix_L, diag_D, matrix_U);
+
+    /*
+        Eliminate lower left triangle.
+    */
+    // For each column,
+    for (int j = 0; j < N - 1; ++j)
+    {
+        const void* block_j = Recovery[j]->Block;
+
+        // For each row,
+        for (int i = j + 1; i < N; ++i)
+        {
+            void* block_i = Recovery[i]->Block;
+            const uint8_t c_ij = *matrix_L++; // Matrix elements are stored column-first, top-down.
+
+            gf256_muladd_mem(block_i, c_ij, block_j, Params.BlockBytes);
+        }
+    }
+
+    /*
+        Eliminate diagonal.
+    */
+    for (int i = 0; i < N; ++i)
+    {
+        void* block = Recovery[i]->Block;
+
+        Recovery[i]->Index = ErasuresIndices[i];
+
+        gf256_div_mem(block, block, diag_D[i], Params.BlockBytes);
+    }
+
+    /*
+        Eliminate upper right triangle.
+    */
+    for (int j = N - 1; j >= 1; --j)
+    {
+        const void* block_j = Recovery[j]->Block;
+
+        for (int i = j - 1; i >= 0; --i)
+        {
+            void* block_i = Recovery[i]->Block;
+            const uint8_t c_ij = *matrix_U++; // Matrix elements are stored column-first, bottom-up.
+
+            gf256_muladd_mem(block_i, c_ij, block_j, Params.BlockBytes);
+        }
+    }
+
+    delete[] dynamicMatrix;
+}
+
+extern "C" int cm256_decode(
+    cm256_encoder_params params, // Encoder params
+    cm256_block* blocks)         // Array of 'originalCount' blocks as described above
+{
+    if (params.OriginalCount <= 0 ||
+        params.RecoveryCount <= 0 ||
+        params.BlockBytes <= 0)
+    {
+        return -1;
+    }
+    if (params.OriginalCount + params.RecoveryCount > 256)
+    {
+        return -2;
+    }
+    if (!blocks)
+    {
+        return -3;
+    }
+
+    // If there is only one block,
+    if (params.OriginalCount == 1)
+    {
+        // It is the same block repeated
+        blocks[0].Index = 0;
+        return 0;
+    }
+
+    CM256Decoder state;
+    if (!state.Initialize(params, blocks))
+    {
+        return -5;
+    }
+
+    // If nothing is erased,
+    if (state.RecoveryCount <= 0)
+    {
+        return 0;
+    }
+
+    // If m=1,
+    if (params.RecoveryCount == 1)
+    {
+        state.DecodeM1();
+        return 0;
+    }
+
+    // Decode for m>1
+    state.Decode();
+    return 0;
+}
diff --git a/gf256.cpp b/gf256.cpp
index 510146b..cad3a43 100644
--- a/gf256.cpp
+++ b/gf256.cpp
@@ -27,7 +27,7 @@
     POSSIBILITY OF SUCH DAMAGE.
 */
 
-#include "gf256.h"
+#include <wirehair/gf256.h>
 
 #ifdef LINUX_ARM
 #include <unistd.h>
diff --git a/include/wirehair/cm256.h b/include/wirehair/cm256.h
new file mode 100644
index 0000000..72a8415
--- /dev/null
+++ b/include/wirehair/cm256.h
@@ -0,0 +1,170 @@
+/*
+	Copyright (c) 2015 Christopher A. Taylor.  All rights reserved.
+
+	Redistribution and use in source and binary forms, with or without
+	modification, are permitted provided that the following conditions are met:
+
+	* Redistributions of source code must retain the above copyright notice,
+	  this list of conditions and the following disclaimer.
+	* Redistributions in binary form must reproduce the above copyright notice,
+	  this list of conditions and the following disclaimer in the documentation
+	  and/or other materials provided with the distribution.
+	* Neither the name of CM256 nor the names of its contributors may be
+	  used to endorse or promote products derived from this software without
+	  specific prior written permission.
+
+	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+	AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+	ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+	LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+	SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+	INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+	CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+	ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+	POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef CM256_H
+#define CM256_H
+
+#include "gf256.h"
+
+#include <assert.h>
+
+// Library version
+#define CM256_VERSION 2
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*
+ * Verify binary compatibility with the API on startup.
+ *
+ * Example:
+ * 	if (cm256_init()) exit(1);
+ *
+ * Returns 0 on success, and any other code indicates failure.
+ */
+extern int cm256_init_(int version);
+#define cm256_init() cm256_init_(CM256_VERSION)
+
+
+// Encoder parameters
+typedef struct cm256_encoder_params_t {
+    // Original block count < 256
+    int OriginalCount;
+
+    // Recovery block count < 256
+    int RecoveryCount;
+
+    // Number of bytes per block (all blocks are the same size in bytes)
+    int BlockBytes;
+} cm256_encoder_params;
+
+// Descriptor for data block
+typedef struct cm256_block_t {
+    // Pointer to data received.
+    void* Block;
+
+    // Block index.
+    // For original data, it will be in the range
+    //    [0..(originalCount-1)] inclusive.
+    // For recovery data, the first one's Index must be originalCount,
+    //    and it will be in the range
+    //    [originalCount..(originalCount+recoveryCount-1)] inclusive.
+    unsigned char Index;
+    // Ignored during encoding, required during decoding.
+} cm256_block;
+
+
+// Compute the value to put in the Index member of cm256_block
+static inline unsigned char cm256_get_recovery_block_index(cm256_encoder_params params, int recoveryBlockIndex)
+{
+    assert(recoveryBlockIndex >= 0 && recoveryBlockIndex < params.RecoveryCount);
+    return (unsigned char)(params.OriginalCount + recoveryBlockIndex);
+}
+static inline unsigned char cm256_get_original_block_index(cm256_encoder_params params, int originalBlockIndex)
+{
+    assert(originalBlockIndex >= 0 && originalBlockIndex < params.OriginalCount);
+    return (unsigned char)(originalBlockIndex);
+}
+
+
+/*
+ * Cauchy MDS GF(256) encode
+ *
+ * This produces a set of recovery blocks that should be transmitted after the
+ * original data blocks.
+ *
+ * It takes in 'originalCount' equal-sized blocks and produces 'recoveryCount'
+ * equally-sized recovery blocks.
+ *
+ * The input 'originals' array allows more natural usage of the library.
+ * The output recovery blocks are stored end-to-end in 'recoveryBlocks'.
+ * 'recoveryBlocks' should have recoveryCount * blockBytes bytes available.
+ *
+ * Precondition: originalCount + recoveryCount <= 256
+ *
+ * When transmitting the data, the block index of the data should be sent,
+ * and the recovery block index is also needed.  The decoder should also
+ * be provided with the values of originalCount, recoveryCount and blockBytes.
+ *
+ * Example wire format:
+ * [originalCount(1 byte)] [recoveryCount(1 byte)]
+ * [blockIndex(1 byte)] [blockData(blockBytes bytes)]
+ *
+ * Be careful not to mix blocks from different encoders.
+ *
+ * It is possible to support variable-length data by including the original
+ * data length at the front of each message in 2 bytes, such that when it is
+ * recovered after a loss the data length is available in the block data and
+ * the remaining bytes of padding can be neglected.
+ *
+ * Returns 0 on success, and any other code indicates failure.
+ */
+extern int cm256_encode(
+    cm256_encoder_params params, // Encoder parameters
+    cm256_block* originals,      // Array of pointers to original blocks
+    void* recoveryBlocks);       // Output recovery blocks end-to-end
+
+// Encode one block.
+// Note: This function does not validate input, use with care.
+extern void cm256_encode_block(
+    cm256_encoder_params params, // Encoder parameters
+    cm256_block* originals,      // Array of pointers to original blocks
+    int recoveryBlockIndex,      // Return value from cm256_get_recovery_block_index()
+    void* recoveryBlock);        // Output recovery block
+
+/*
+ * Cauchy MDS GF(256) decode
+ *
+ * This recovers the original data from the recovery data in the provided
+ * blocks.  There should be 'originalCount' blocks in the provided array.
+ * Recovery will always be possible if that many blocks are received.
+ *
+ * Provide the same values for 'originalCount', 'recoveryCount', and
+ * 'blockBytes' used by the encoder.
+ *
+ * The block Index should be set to the block index of the original data,
+ * as described in the cm256_block struct comments above.
+ *
+ * Recovery blocks will be replaced with original data and the Index
+ * will be updated to indicate the original block that was recovered.
+ *
+ * Returns 0 on success, and any other code indicates failure.
+ */
+extern int cm256_decode(
+    cm256_encoder_params params, // Encoder parameters
+    cm256_block* blocks);        // Array of 'originalCount' blocks as described above
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif // CM256_H
diff --git a/gf256.h b/include/wirehair/gf256.h
similarity index 100%
rename from gf256.h
rename to include/wirehair/gf256.h
diff --git a/tables/HeavyRowGenerator.cpp b/tables/HeavyRowGenerator.cpp
index c8ec0c8..c3050a6 100644
--- a/tables/HeavyRowGenerator.cpp
+++ b/tables/HeavyRowGenerator.cpp
@@ -36,7 +36,7 @@ using namespace std;
 static const unsigned kRows = 6;
 static const unsigned kColumns = 18;
 
-#include "../gf256.h"
+#include <wirehair/gf256.h>
 #include "../test/SiameseTools.h"
 using namespace siamese;
 
