An implementation of the finite field of order 2^256.

f_2_256.cc

/*

An implementation of the finite field of order 2^256.

The MIT License (MIT)
Copyright (c) 2024 Ji Luo

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

cl /EHsc /std:c++17 f_2_256.cc

*/

#include<utility>
#include<cstring>

/*
The irreducible polynomial is
  p(x) = x^(256-0) + x^(256-15) + x^(256-78) + x^(256-135) + x^(256-256).
*/
struct GF_2_256
{
  unsigned char data[32];
  constexpr GF_2_256() noexcept = default;
  constexpr GF_2_256(GF_2_256 const &) noexcept = default;
  constexpr GF_2_256(GF_2_256 &&) noexcept = default;
  explicit GF_2_256(unsigned char const *data_) noexcept
  {
    std::memcpy(data, data_, 32);
  }
  constexpr GF_2_256 &operator =(GF_2_256 const &) noexcept = default;
  constexpr GF_2_256 &operator =(GF_2_256 &&) noexcept = default;
  GF_2_256 &operator =(unsigned char const *data_) noexcept
  {
    std::memmove(data, data_, 32);
    return *this;
  }
  ~GF_2_256() noexcept = default;
private:
  struct tag1 { };
  constexpr GF_2_256(tag1) noexcept : data{1}
  {
  }
public:
  static constexpr GF_2_256 zero() { return {}; }
  static constexpr GF_2_256 one() { return tag1{}; }
  /* positive */
  friend GF_2_256 const &operator +(GF_2_256 const &a)
  {
    return a;
  }
  friend GF_2_256 &&operator +(GF_2_256 &&a)
  {
    return std::move(a);
  }
  friend GF_2_256 const &operator -(GF_2_256 const &a)
  {
    return a;
  }
  friend GF_2_256 &&operator -(GF_2_256 &&a)
  {
    return std::move(a);
  }
  /* addition */
  friend GF_2_256 operator +(GF_2_256 const &a, GF_2_256 const &b)
  {
    GF_2_256 result{};
    for (size_t i = 0; i != 32; ++i)
      result.data[i] = a.data[i] ^ b.data[i];
    return result;
  }
  friend GF_2_256 operator -(GF_2_256 const &a, GF_2_256 const &b)
  {
    return a + b;
  }
  /* multiplication */
private:
  struct impl
  {
    /*
    inline static constexpr unsigned char product[256][256][2] = { ... };
    inline static constexpr unsigned char remainder[32][256][32] = { ... };

    product[r][s] = r(x) s(x)
    remainder[t][q] = x^(8t) q(x) mod p(x)
    */
    #include"./f_2_256_computed.cc"
  };
public:
  friend GF_2_256 operator *(GF_2_256 const &a, GF_2_256 const &b)
  {
    unsigned char result[64];
    for (size_t i = 0; i != 32; ++i)
      for (size_t j = 0; j != 32; ++j)
        result[i + j] ^= impl::product[a.data[i]][b.data[j]][0],
        result[i + j + 1] ^= impl::product[a.data[i]][b.data[j]][1];
    for (size_t t = 0; t != 32; ++t)
      for (size_t k = 0; k != 32; ++k)
        result[k] ^= impl::remainder[t][result[32 + t]][k];
    return GF_2_256{result};
  }
  /* division */
  static GF_2_256 reciprocal(GF_2_256 a)
  {
    /* This is a stupid algorithm... */
    GF_2_256 result{tag1{}};
    for (size_t i = 1; i != 256; ++i)
    {
      a = a * a;
      result = result * a;
    }
    return result;
  }
  GF_2_256 reciprocal() const
  {
    return reciprocal(*this);
  }
  friend GF_2_256 operator /(GF_2_256 const &a, GF_2_256 b)
  {
    return a * reciprocal(b);
  }
  /* compound assignment operators */
  GF_2_256 &operator +=(GF_2_256 const &that)
  {
    for (size_t i = 0; i != 32; ++i)
      data[i] ^= that.data[i];
    return *this;
  }
  GF_2_256 &operator -=(GF_2_256 const &that)
  {
    return *this += that;
  }
  GF_2_256 &operator *=(GF_2_256 const &that)
  {
    return *this = *this * that;
  }
  GF_2_256 &operator /=(GF_2_256 that)
  {
    return *this = *this * reciprocal(that);
  }
};

int main()
{
}

f_2_256_compute.cc

/*

An implementation of the finite field of order 2^256.

The MIT License (MIT)
Copyright (c) 2024 Ji Luo

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

cl /EHsc /std:c++17 f_2_256_compute.cc && f_2_256_compute

*/

#include<vector>
#include<utility>

struct F_2_x
{
  std::vector<bool> coefficients;
  F_2_x() = default;
  F_2_x(F_2_x const &) = default;
  F_2_x(F_2_x &&) = default;
  explicit F_2_x(std::vector<bool> coefficients_)
    : coefficients{ std::move(coefficients_) }
  {
  }
  explicit F_2_x(unsigned char coefficients_)
  {
    size_t s = (coefficients_ > 127u ? (size_t)8
      : coefficients_ > 63u ? (size_t)7
      : coefficients_ > 31u ? (size_t)6
      : coefficients_ > 15u ? (size_t)5
      : coefficients_ > 7u ? (size_t)4
      : coefficients_ > 3u ? (size_t)3
      : coefficients_ > 1u ? (size_t)2
      : coefficients_ > 0u ? (size_t)1
      : (size_t)0);
    coefficients.resize(s);
    for (size_t i = (size_t)0; i != s; ++i)
      coefficients[i] = (coefficients_ & (1u << i));
  }
  F_2_x &operator =(F_2_x const &) = default;
  F_2_x &operator =(F_2_x &&) = default;
  F_2_x &operator =(std::vector<bool> const &coefficients_)
  {
    coefficients = coefficients_;
    return *this;
  }
  F_2_x &operator =(std::vector<bool> &&coefficients_)
  {
    coefficients = std::move(coefficients_);
    return *this;
  }
  ~F_2_x() = default;
  friend F_2_x operator +(F_2_x a)
  {
    size_t i = a.coefficients.size();
    while (i && !a.coefficients[i - (size_t)1])
      --i;
    a.coefficients.resize(i);
    return a;
  }
  friend F_2_x operator +(F_2_x const &a, F_2_x const &b)
  {
    F_2_x result;
    size_t aa = a.coefficients.size(), bb = b.coefficients.size();
    result.coefficients.resize(std::max(aa, bb));
    size_t i = (size_t)0;
    for (size_t mm = std::min(aa, bb); i != mm; ++i)
      result.coefficients[i] = a.coefficients[i] ^ b.coefficients[i];
    if (aa > bb)
      for (; i != aa; ++i)
        result.coefficients[i] = a.coefficients[i];
    if (aa < bb)
      for (; i != bb; ++i)
        result.coefficients[i] = b.coefficients[i];
    while (i && !result.coefficients[i - (size_t)1])
      --i;
    result.coefficients.resize(i);
    return result;
  }
  friend F_2_x operator *(F_2_x const &a, F_2_x const &b)
  {
    size_t aa = a.coefficients.size();
    while (aa && !a.coefficients[aa - (size_t)1])
      --aa;
    size_t bb = b.coefficients.size();
    while (bb && !b.coefficients[bb - (size_t)1])
      --bb;
    F_2_x result;
    if (aa && bb)
    {
      result.coefficients.resize(aa + bb - (size_t)1);
      for (size_t i = (size_t)0; i != aa; ++i)
        if (a.coefficients[i])
          for (size_t j = (size_t)0; j != bb; ++j)
            result.coefficients[i + j] = result.coefficients[i + j] ^ b.coefficients[j];
    }
    return result;
  }
  friend F_2_x operator %(F_2_x a, F_2_x const &b)
  {
    size_t aa = a.coefficients.size();
    while (aa && !a.coefficients[aa - (size_t)1])
      --aa;
    size_t bb = b.coefficients.size();
    while (bb && !b.coefficients[bb - (size_t)1])
      --bb;
    /* Division by zero is supposed to be undefined,
    /* which happens to be implemented by an infinite loop,
    /* which by the C++ standard is also undefined. */
    while (aa >= bb)
    {
      for (size_t i = (size_t)0, d = aa - bb; i != bb; ++i)
        a.coefficients[i + d] = a.coefficients[i + d] ^ b.coefficients[i];
      while (aa && !a.coefficients[aa - (size_t)1])
        --aa;
    }
    a.coefficients.resize(aa);
    return a;
  }
  unsigned char get(size_t from) const
  {
    unsigned char result{};
    size_t to = std::min(from + (size_t)8, coefficients.size());
    for (size_t i = from; i < to; ++i)
      result |= (coefficients[i] << (i - from));
    return result;
  }
};

#include<cstdio>

void compute_products()
{
  std::fputs("inline static constexpr unsigned char product[256][256][2] = {",
    stdout);
  bool first = true;
  for (unsigned r = 0u; r != 256u; ++r)
  {
    std::fputs(first ? " {" : ", {", stdout);
    first = true;
    for (unsigned s = 0u; s != 256u; ++s)
    {
      F_2_x product = F_2_x{(unsigned char)r} * F_2_x{(unsigned char)s};
      std::printf("%s{ %u, %u }",
        first ? " " : ", ",
        (unsigned)product.get((size_t)0),
        (unsigned)product.get((size_t)8));
      first = false;
    }
    std::fputs(" }", stdout);
    first = false;
  }
  std::fputs(" };\n", stdout);
}

void compute_remainders()
{
  std::fputs("inline static constexpr unsigned char remainder[32][256][32] = {", stdout);
  bool first = true;
  F_2_x p;
  /*
  The irreducible polynomial is
    p(x) = x^(256-0) + x^(256-15) + x^(256-78) + x^(256-135) + x^(256-256).
  */
  p.coefficients.resize(257);
  p.coefficients[256 - 0] = true;
  p.coefficients[256 - 15] = true;
  p.coefficients[256 - 78] = true;
  p.coefficients[256 - 135] = true;
  p.coefficients[256 - 256] = true;
  for (unsigned t = 0u; t != 32u; ++t)
  {
    std::fputs(first ? " {" : ", {", stdout);
    first = true;
    F_2_x x8t;
    x8t.coefficients.resize(8 * t + 1);
    x8t.coefficients[8 * t] = true;
    for (unsigned q = 0u; q != 256u; ++q)
    {
      std::fputs(first ? " {" : ", {", stdout);
      first = true;
      F_2_x remainder = (x8t * F_2_x{(unsigned char)q}) % p;
      for (unsigned part = 0u; part != 32u; ++part)
      {
        std::printf("%s%u",
          first ? " " : ", ",
          (unsigned)remainder.get((size_t)(part * 8u)));
        first = false;
      }
      std::fputs(" }", stdout);
      first = false;
    }
    std::fputs(" }", stdout);
    first = false;
  }
  std::fputs(" };\n", stdout);
}

int main()
{
  std::freopen("f_2_256_computed.cc", "wb", stdout);
  compute_products();
  compute_remainders();
}