Skip to content

description: Compute the 64-bit nim product. Nimbers with twice as many bits form a quadratic extension of the lower field: for alpha=2^m, alpha^2=alpha xor 2^(m-1). Recursing on this identity generates all products of binary basis elements. Bilinearity then groups eight basis elements at a time into a 32 MiB byte-pair table, so one query needs 8*8 lookups and XORs.

nim_product.hpp

SECTIONMath INCLUDEnoya/nim_product.hpp

Compute the 64-bit nim product. Nimbers with twice as many bits form a quadratic extension of the lower field: for alpha=2^m, alpha^2=alpha xor 2^(m-1). Recursing on this identity generates all products of binary basis elements. Bilinearity then groups eight basis elements at a time into a 32 MiB byte-pair table, so one query needs 8*8 lookups and XORs.

Verified by nim_product_64.

\[ \displaystyle a \otimes b \]

Implementation

View on GitHub

#ifndef NOYA_NIM_PRODUCT_HPP
#define NOYA_NIM_PRODUCT_HPP 1

/// @complexity Time: O(64) per product after O(2^22) preprocessing.
/// Space: O(2^25) bytes (32 MiB).

#include <array>
#include <bit>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <vector>

namespace noya {

namespace nim_product_internal {

using u64 = std::uint64_t;

inline u64 recursive_product(u64 first, u64 second, int bit_width) {
  assert(bit_width >= 1 && bit_width <= 64 &&
         (bit_width & (bit_width - 1)) == 0);
  if (first == 0 || second == 0) {
    return 0;
  }
  if (first == 1) {
    return second;
  }
  if (second == 1) {
    return first;
  }
  if (bit_width == 1) {
    return first & second;
  }

  int half = bit_width / 2;
  u64 mask = (u64(1) << half) - 1;
  u64 first_low = first & mask;
  u64 first_high = first >> half;
  u64 second_low = second & mask;
  u64 second_high = second >> half;

  u64 low_low = recursive_product(first_low, second_low, half);
  u64 high_high = recursive_product(first_high, second_high, half);
  u64 combined = recursive_product(first_low ^ first_high,
                                   second_low ^ second_high, half);
  u64 cross = combined ^ low_low ^ high_high;
  u64 reduction = recursive_product(high_high, u64(1) << (half - 1), half);
  return (low_low ^ reduction) ^ ((cross ^ high_high) << half);
}

inline const std::vector<u64> &byte_table() {
  static const std::vector<u64> table = [] {
    std::array<std::array<u64, 64>, 64> basis{};
    for (int first = 0; first < 64; first++) {
      for (int second = first; second < 64; second++) {
        u64 value = recursive_product(u64(1) << first,
                                      u64(1) << second, 64);
        basis[first][second] = basis[second][first] = value;
      }
    }

    std::vector<u64> result(std::size_t(8) * 8 * 256 * 256);
    std::array<std::array<u64, 8>, 256> row_contribution{};
    for (int first_byte = 0; first_byte < 8; first_byte++) {
      for (int second_byte = 0; second_byte < 8; second_byte++) {
        row_contribution = {};
        for (unsigned mask = 1; mask < 256; mask++) {
          unsigned bit = std::countr_zero(mask);
          unsigned previous = mask & (mask - 1);
          for (int second_bit = 0; second_bit < 8; second_bit++) {
            row_contribution[mask][second_bit] =
                row_contribution[previous][second_bit] ^
                basis[first_byte * 8 + int(bit)]
                     [second_byte * 8 + second_bit];
          }
        }
        std::size_t block =
            std::size_t(first_byte * 8 + second_byte) << 16;
        for (unsigned first_value = 0; first_value < 256; first_value++) {
          for (unsigned second_value = 1; second_value < 256;
               second_value++) {
            unsigned bit = std::countr_zero(second_value);
            unsigned previous = second_value & (second_value - 1);
            result[block | (std::size_t(first_value) << 8) | second_value] =
                result[block | (std::size_t(first_value) << 8) | previous] ^
                row_contribution[first_value][bit];
          }
        }
      }
    }
    return result;
  }();
  return table;
}

} // namespace nim_product_internal

/// @brief Compute the 64-bit nim product. Nimbers with twice as many bits form
/// a quadratic extension of the lower field: for alpha=2^m,
/// alpha^2=alpha xor 2^(m-1). Recursing on this identity generates all products
/// of binary basis elements. Bilinearity then groups eight basis elements at a
/// time into a 32 MiB byte-pair table, so one query needs 8*8 lookups and XORs.
inline std::uint64_t nim_product(std::uint64_t first, std::uint64_t second) {
  const std::vector<std::uint64_t> &table =
      nim_product_internal::byte_table();
  std::uint64_t result = 0;
  for (int first_byte = 0; first_byte < 8; first_byte++) {
    unsigned first_value = unsigned(first >> (first_byte * 8)) & 255;
    for (int second_byte = 0; second_byte < 8; second_byte++) {
      unsigned second_value = unsigned(second >> (second_byte * 8)) & 255;
      std::size_t index =
          (std::size_t(first_byte * 8 + second_byte) << 16) |
          (std::size_t(first_value) << 8) | second_value;
      result ^= table[index];
    }
  }
  return result;
}

} // namespace noya

#endif // NOYA_NIM_PRODUCT_HPP
#include <array>
#include <bit>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <vector>

/// @complexity Time: O(64) per product after O(2^22) preprocessing.
/// Space: O(2^25) bytes (32 MiB).

namespace noya {

namespace nim_product_internal {

using u64 = std::uint64_t;

inline u64 recursive_product(u64 first, u64 second, int bit_width) {
  assert(bit_width >= 1 && bit_width <= 64 &&
         (bit_width & (bit_width - 1)) == 0);
  if (first == 0 || second == 0) {
    return 0;
  }
  if (first == 1) {
    return second;
  }
  if (second == 1) {
    return first;
  }
  if (bit_width == 1) {
    return first & second;
  }

  int half = bit_width / 2;
  u64 mask = (u64(1) << half) - 1;
  u64 first_low = first & mask;
  u64 first_high = first >> half;
  u64 second_low = second & mask;
  u64 second_high = second >> half;

  u64 low_low = recursive_product(first_low, second_low, half);
  u64 high_high = recursive_product(first_high, second_high, half);
  u64 combined = recursive_product(first_low ^ first_high,
                                   second_low ^ second_high, half);
  u64 cross = combined ^ low_low ^ high_high;
  u64 reduction = recursive_product(high_high, u64(1) << (half - 1), half);
  return (low_low ^ reduction) ^ ((cross ^ high_high) << half);
}

inline const std::vector<u64> &byte_table() {
  static const std::vector<u64> table = [] {
    std::array<std::array<u64, 64>, 64> basis{};
    for (int first = 0; first < 64; first++) {
      for (int second = first; second < 64; second++) {
        u64 value = recursive_product(u64(1) << first,
                                      u64(1) << second, 64);
        basis[first][second] = basis[second][first] = value;
      }
    }

    std::vector<u64> result(std::size_t(8) * 8 * 256 * 256);
    std::array<std::array<u64, 8>, 256> row_contribution{};
    for (int first_byte = 0; first_byte < 8; first_byte++) {
      for (int second_byte = 0; second_byte < 8; second_byte++) {
        row_contribution = {};
        for (unsigned mask = 1; mask < 256; mask++) {
          unsigned bit = std::countr_zero(mask);
          unsigned previous = mask & (mask - 1);
          for (int second_bit = 0; second_bit < 8; second_bit++) {
            row_contribution[mask][second_bit] =
                row_contribution[previous][second_bit] ^
                basis[first_byte * 8 + int(bit)]
                     [second_byte * 8 + second_bit];
          }
        }
        std::size_t block =
            std::size_t(first_byte * 8 + second_byte) << 16;
        for (unsigned first_value = 0; first_value < 256; first_value++) {
          for (unsigned second_value = 1; second_value < 256;
               second_value++) {
            unsigned bit = std::countr_zero(second_value);
            unsigned previous = second_value & (second_value - 1);
            result[block | (std::size_t(first_value) << 8) | second_value] =
                result[block | (std::size_t(first_value) << 8) | previous] ^
                row_contribution[first_value][bit];
          }
        }
      }
    }
    return result;
  }();
  return table;
}

} // namespace nim_product_internal

/// @brief Compute the 64-bit nim product. Nimbers with twice as many bits form
/// a quadratic extension of the lower field: for alpha=2^m,
/// alpha^2=alpha xor 2^(m-1). Recursing on this identity generates all products
/// of binary basis elements. Bilinearity then groups eight basis elements at a
/// time into a 32 MiB byte-pair table, so one query needs 8*8 lookups and XORs.
inline std::uint64_t nim_product(std::uint64_t first, std::uint64_t second) {
  const std::vector<std::uint64_t> &table =
      nim_product_internal::byte_table();
  std::uint64_t result = 0;
  for (int first_byte = 0; first_byte < 8; first_byte++) {
    unsigned first_value = unsigned(first >> (first_byte * 8)) & 255;
    for (int second_byte = 0; second_byte < 8; second_byte++) {
      unsigned second_value = unsigned(second >> (second_byte * 8)) & 255;
      std::size_t index =
          (std::size_t(first_byte * 8 + second_byte) << 16) |
          (std::size_t(first_value) << 8) | second_value;
      result ^= table[index];
    }
  }
  return result;
}

} // namespace noya