symmetryimpl.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <uni20/common/trace.hpp>
#include <uni20/symmetry/u1.hpp>
 
#include <cctype>
#include <cstdint>
#include <mutex>
#include <optional>
#include <span>
#include <sstream>
#include <stdexcept>
#include <string>
#include <string_view>
#include <unordered_map>
#include <utility>
#include <vector>
 
namespace uni20::detail
{
 
inline constexpr std::uint64_t qnum_out_of_line_mask = std::uint64_t{1} << 63;
 
inline constexpr std::uint64_t qnum_inline_code_max = qnum_out_of_line_mask - 1;
 
struct SymmetryFactorSpec
{
    std::string name;
    SymmetryFactorBase const* factor = nullptr;
};
 
class SymmetryImpl {
  public:
    explicit SymmetryImpl(std::vector<SymmetryFactorSpec> factors) : factors_(std::move(factors))
    {
      for (std::size_t i = 0; i < factors_.size(); ++i)
      {
        auto const [it, inserted] = index_by_name_.emplace(factors_[i].name, i);
        if (!inserted)
        {
          throw std::invalid_argument("duplicate symmetry component name: " + it->first);
        }
      }
      canonical_string_ = build_canonical_string(factors_);
    }
 
    auto factors() const -> std::span<SymmetryFactorSpec const> { return factors_; }
 
    auto factor_count() const -> std::size_t { return factors_.size(); }
 
    auto find_factor(std::string_view name) const -> std::optional<std::size_t>
    {
      auto const it = index_by_name_.find(std::string{name});
      if (it == index_by_name_.end())
      {
        return std::nullopt;
      }
      return it->second;
    }
 
    auto to_string() const -> std::string const& { return canonical_string_; }
 
    static auto intern(std::vector<SymmetryFactorSpec> factors) -> SymmetryImpl const*
    {
      auto const key = build_canonical_string(factors);
 
      auto& mutex = intern_mutex();
      std::lock_guard<std::mutex> lock(mutex);
 
      auto& map = intern_map();
      if (auto const it = map.find(key); it != map.end())
      {
        return it->second;
      }
 
      auto* impl = new SymmetryImpl(std::move(factors));
      map.emplace(impl->to_string(), impl);
      return impl;
    }
 
    static auto parse(std::string_view spec) -> SymmetryImpl const*
    {
      register_builtin_symmetry_factors();
 
      auto trimmed = trim(spec);
      if (trimmed.empty())
      {
        throw std::invalid_argument("symmetry specification must not be empty");
      }
 
      std::vector<SymmetryFactorSpec> factors;
      std::size_t cursor = 0;
      while (cursor < trimmed.size())
      {
        auto const next = trimmed.find(',', cursor);
        auto const token =
            trim(trimmed.substr(cursor, next == std::string_view::npos ? std::string_view::npos : next - cursor));
        if (token.empty())
        {
          throw std::invalid_argument("symmetry specification contains an empty factor");
        }
 
        auto const colon = token.find(':');
        if (colon == std::string_view::npos)
        {
          throw std::invalid_argument("symmetry factor is missing ':' separator: " + std::string{token});
        }
 
        auto const name = trim(token.substr(0, colon));
        auto const kind = trim(token.substr(colon + 1));
        if (name.empty())
        {
          throw std::invalid_argument("symmetry factor has an empty name");
        }
        if (kind.empty())
        {
          throw std::invalid_argument("symmetry factor has an empty type");
        }
 
        factors.push_back({std::string{name}, find_symmetry_factor(kind)});
 
        if (next == std::string_view::npos)
        {
          break;
        }
        cursor = next + 1;
      }
 
      return intern(std::move(factors));
    }
 
    auto pack(std::span<std::uint64_t const> factor_codes) const -> std::uint64_t
    {
      if (factor_codes.size() != factors_.size())
      {
        throw std::invalid_argument("wrong number of factor codes for symmetry packing");
      }
      if (factor_codes.empty())
      {
        return 0;
      }
 
      std::uint64_t combined = factor_codes.front();
      if (combined > qnum_inline_code_max)
      {
        throw std::overflow_error("quantum number code requires the reserved out-of-line representation");
      }
      for (std::size_t i = 1; i < factor_codes.size(); ++i)
      {
        combined = cantor_pair(combined, factor_codes[i]);
      }
      return combined;
    }
 
    auto unpack(std::uint64_t code) const -> std::vector<std::uint64_t>
    {
      if (code > qnum_inline_code_max)
      {
        throw std::invalid_argument("cannot unpack a tagged out-of-line quantum number as an inline code");
      }
      if (factors_.empty())
      {
        return {};
      }
      if (factors_.size() == 1)
      {
        return {code};
      }
 
      std::vector<std::uint64_t> codes(factors_.size());
      auto current = code;
      for (std::size_t i = factors_.size(); i-- > 1;)
      {
        auto const [left, right] = cantor_unpair(current);
        codes[i] = right;
        current = left;
      }
      codes[0] = current;
      return codes;
    }
 
  private:
    static auto trim(std::string_view text) -> std::string_view
    {
      while (!text.empty() && std::isspace(static_cast<unsigned char>(text.front())) != 0)
      {
        text.remove_prefix(1);
      }
      while (!text.empty() && std::isspace(static_cast<unsigned char>(text.back())) != 0)
      {
        text.remove_suffix(1);
      }
      return text;
    }
 
    static auto build_canonical_string(std::span<SymmetryFactorSpec const> factors) -> std::string
    {
      std::ostringstream out;
      for (std::size_t i = 0; i < factors.size(); ++i)
      {
        if (i != 0)
        {
          out << ',';
        }
        out << factors[i].name << ':' << factors[i].factor->type_name();
      }
      return out.str();
    }
 
    static auto intern_map() -> std::unordered_map<std::string, SymmetryImpl const*>&
    {
      // Canonical instances have process lifetime. Keep their registry alive
      // too, so static destruction does not orphan the retained instances.
      static auto* const map = new std::unordered_map<std::string, SymmetryImpl const*>;
      return *map;
    }
 
    static auto intern_mutex() -> std::mutex&
    {
      static std::mutex mutex;
      return mutex;
    }
 
    static auto triangular_leq(std::uint64_t value, std::uint64_t limit) -> bool
    {
      auto const left = (value % 2U == 0U) ? value / 2 : value / 2 + 1;
      auto const right = (value % 2U == 0U) ? value + 1 : value;
      if (left == 0)
      {
        return true;
      }
      return left <= limit / right;
    }
 
    static auto triangular(std::uint64_t value) -> std::uint64_t
    {
      auto const left = (value % 2U == 0U) ? value / 2 : value / 2 + 1;
      auto const right = (value % 2U == 0U) ? value + 1 : value;
      return left * right;
    }
 
    static auto cantor_pair(std::uint64_t left, std::uint64_t right) -> std::uint64_t
    {
      if (left > qnum_inline_code_max || right > qnum_inline_code_max || left > qnum_inline_code_max - right)
      {
        throw std::overflow_error("quantum number code requires the reserved out-of-line representation");
      }
 
      auto const sum = left + right;
      if (!triangular_leq(sum, qnum_inline_code_max - right))
      {
        throw std::overflow_error("quantum number code requires the reserved out-of-line representation");
      }
      return triangular(sum) + right;
    }
 
    static auto cantor_unpair(std::uint64_t code) -> std::pair<std::uint64_t, std::uint64_t>
    {
      std::uint64_t lo = 0;
      std::uint64_t hi = code;
      while (lo < hi)
      {
        auto const mid = lo + (hi - lo + 1) / 2;
        if (triangular_leq(mid, code))
        {
          lo = mid;
        }
        else
        {
          hi = mid - 1;
        }
      }
 
      auto const sum = lo;
      auto const triangle = triangular(sum);
      auto const right = code - triangle;
      auto const left = sum - right;
      return {left, right};
    }
 
    std::vector<SymmetryFactorSpec> factors_;
    std::unordered_map<std::string, std::size_t> index_by_name_;
    std::string canonical_string_;
};
 
} // namespace uni20::detail