Suffix_Structures.H

Go to the documentation of this file.

/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  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.
*/
 
#ifndef SUFFIX_STRUCTURES_H
#define SUFFIX_STRUCTURES_H
 
#include <algorithm>
#include <limits>
#include <string>
#include <string_view>
#include <utility>
#include <array>
 
#include <ah-errors.H>
#include <tpl_array.H>
#include <tpl_sort_utils.H>
 
namespace Aleph {
namespace suffix_structures_detail {
inline Array<size_t> all_positions(const size_t n)
{
  Array<size_t> out;
  out.reserve(n + 1);
  for (size_t i = 0; i <= n; ++i)
    out.append(i);
  return out;
}
 
inline void sort_matches(Array<size_t> &a)
{
  introsort(a);
}
}  // namespace suffix_structures_detail
 
inline Array<size_t> suffix_array(const std::string_view text)
{
  const size_t n = text.size();
  if (n == 0)
    return {};
 
  Array<size_t> sa = Array<size_t>::create(n);
  Array<int> rank = Array<int>::create(n);
  Array<int> tmp_rank = Array<int>::create(n);
  Array<size_t> tmp_sa = Array<size_t>::create(n);
 
  for (size_t i = 0; i < n; ++i)
    {
      sa[i] = i;
      rank[i] = static_cast<unsigned char>(text[i]);
    }
 
  int max_rank = 255;  // initial ranks are char codes in [0, 255]
 
  for (size_t k = 1;;)
    {
      auto key2 = [&](const size_t idx) -> int
      {
        return (idx + k < n) ? rank[idx + k] : -1;
      };
 
      auto key1 = [&](const size_t idx) -> int
      {
        return rank[idx];
      };
 
      counting_sort_indices(sa, tmp_sa, n, -1, max_rank, key2);
      counting_sort_indices(sa, tmp_sa, n, 0, max_rank, key1);
 
      auto less_idx = [&](const size_t a, const size_t b)
      {
        if (rank[a] != rank[b])
          return rank[a] < rank[b];
        const int ra = (a + k < n) ? rank[a + k] : -1;
        const int rb = (b + k < n) ? rank[b + k] : -1;
        return ra < rb;
      };
 
      tmp_rank[sa[0]] = 0;
      for (size_t i = 1; i < n; ++i)
        tmp_rank[sa[i]] = tmp_rank[sa[i - 1]] + (less_idx(sa[i - 1], sa[i]) ? 1 : 0);
 
      for (size_t i = 0; i < n; ++i)
        rank[i] = tmp_rank[i];
 
      max_rank = rank[sa[n - 1]];
 
      if (max_rank == static_cast<int>(n - 1))
        break;
 
      if (k >= n)
        break;
 
      if (k > n / 2)
        break;
 
      k *= 2;
    }
 
  return sa;
}
 
inline Array<size_t> lcp_array_kasai(const std::string_view text, const Array<size_t> &sa)
{
  const size_t n = text.size();
 
  ah_domain_error_if(sa.size() != n) << "lcp_array_kasai(): suffix array size must match text size";
 
  Array<size_t> lcp;
  if (n == 0)
    return lcp;
 
  lcp = Array<size_t>::create(n);
  for (size_t i = 0; i < n; ++i)
    lcp[i] = 0;
 
  Array<size_t> rank(n, std::numeric_limits<size_t>::max());
  for (size_t i = 0; i < n; ++i)
    {
      ah_out_of_range_error_if(sa[i] >= n)
        << "lcp_array_kasai(): suffix array contains invalid index";
      ah_out_of_range_error_if(rank[sa[i]] != std::numeric_limits<size_t>::max())
        << "lcp_array_kasai(): suffix array contains duplicate index";
      rank[sa[i]] = i;
    }
 
  size_t k = 0;
  for (size_t i = 0; i < n; ++i)
    {
      const size_t r = rank[i];
      if (r + 1 >= n)
        {
          k = 0;
          continue;
        }
 
      const size_t j = sa[r + 1];
      while (i + k < n and j + k < n and text[i + k] == text[j + k])
        ++k;
 
      lcp[r] = k;
 
      if (k > 0)
        --k;
    }
 
  return lcp;
}
 
class Naive_Suffix_Tree
{
public:
  static constexpr size_t npos = std::numeric_limits<size_t>::max();
 
  struct Node
  {
    size_t start = 0;           
    size_t end = 0;             
    size_t parent = npos;       
    Array<size_t> children;     
    size_t suffix_index = npos; 
  };
 
private:
  std::string text_;
  size_t original_size_ = 0;
  Array<Node> nodes_;
 
  static char choose_terminal(const std::string_view s)
  {
    std::array<bool, 256> used;
    used.fill(false);
    for (const unsigned char c : s)
      used[c] = true;
 
    for (size_t i = 0; i < 256; ++i)
      if (not used[i])
        return static_cast<char>(i);
 
    ah_domain_error() << "Naive_Suffix_Tree::choose_terminal(): all 256 byte values "
                      << "are present in the input; no terminal marker available";
 
    return '\0';  // unreachable
  }
 
  size_t create_node(const size_t start, const size_t end, const size_t parent, const size_t suffix_index)
  {
    Node n;
    n.start = start;
    n.end = end;
    n.parent = parent;
    n.suffix_index = suffix_index;
    nodes_.append(n);
    return nodes_.size() - 1;
  }
 
  size_t find_child_by_first_char(const size_t node, const char c) const
  {
    const Array<size_t> &children = nodes_[node].children;
    for (size_t i = 0; i < children.size(); ++i)
      if (const size_t child = children[i]; text_[nodes_[child].start] == c)
        return child;
 
    return npos;
  }
 
  void replace_child(const size_t parent, const size_t old_child, const size_t new_child)
  {
    for (Array<size_t> &children = nodes_[parent].children; size_t & i : children)
      if (i == old_child)
        {
          i = new_child;
          return;
        }
  }
 
  void insert_suffix(const size_t suffix_start)
  {
    size_t current = 0;
    size_t pos = suffix_start;
 
    while (pos < text_.size())
      {
        const size_t child = find_child_by_first_char(current, text_[pos]);
        if (child == npos)
          {
            const size_t leaf = create_node(pos, text_.size(), current, suffix_start);
            nodes_[current].children.append(leaf);
            return;
          }
 
        const size_t edge_start = nodes_[child].start;
        const size_t edge_end = nodes_[child].end;
 
        size_t k = 0;
        while (pos + k < text_.size() and edge_start + k < edge_end
               and text_[pos + k] == text_[edge_start + k])
          ++k;
 
        if (edge_start + k == edge_end)
          {
            current = child;
            pos += k;
            continue;
          }
 
        const size_t split = create_node(edge_start, edge_start + k, current, npos);
 
        replace_child(current, child, split);
 
        nodes_[child].start = edge_start + k;
        nodes_[child].parent = split;
        nodes_[split].children.append(child);
 
        const size_t leaf = create_node(pos + k, text_.size(), split, suffix_start);
        nodes_[split].children.append(leaf);
 
        return;
      }
  }
 
  void collect_leaf_suffixes(const size_t node, Array<size_t> &out) const
  {
    Array<size_t> stack;
    stack.append(node);
 
    while (not stack.is_empty())
      {
        const size_t u = stack.remove_last();
        if (nodes_[u].suffix_index != npos)
          {
            if (nodes_[u].suffix_index < original_size_)
              out.append(nodes_[u].suffix_index);
            continue;
          }
 
        for (const Array<size_t> &children = nodes_[u].children; size_t i : children)
          stack.append(i);
      }
  }
 
public:
  explicit Naive_Suffix_Tree(const std::string_view text = "")
  {
    build(text);
  }
 
  void build(const std::string_view text)
  {
    const char terminal = choose_terminal(text);
    std::string new_text(text.begin(), text.end());
    new_text.push_back(terminal);
 
    original_size_ = text.size();
    text_ = std::move(new_text);
 
    nodes_.empty();
    nodes_.append(Node{});  // root
 
    for (size_t i = 0; i < text_.size(); ++i)
      insert_suffix(i);
  }
 
  [[nodiscard]] bool contains(const std::string_view pattern) const
  {
    if (pattern.empty())
      return true;
 
    size_t current = 0;
    size_t pos = 0;
 
    while (pos < pattern.size())
      {
        const size_t child = find_child_by_first_char(current, pattern[pos]);
        if (child == npos)
          return false;
 
        const size_t edge_start = nodes_[child].start;
        const size_t edge_end = nodes_[child].end;
 
        size_t k = 0;
        while (pos + k < pattern.size() and edge_start + k < edge_end
               and text_[edge_start + k] != text_.back()  // Never match sentinel
               and pattern[pos + k] == text_[edge_start + k])
          ++k;
 
        if (pos + k == pattern.size())
          return true;
 
        if (edge_start + k == edge_end)
          {
            current = child;
            pos += k;
            continue;
          }
 
        return false;
      }
 
    return true;
  }
 
  [[nodiscard]] Array<size_t> find_all(const std::string_view pattern) const
  {
    if (pattern.empty())
      return suffix_structures_detail::all_positions(original_size_);
 
    size_t current = 0;
    size_t pos = 0;
 
    while (pos < pattern.size())
      {
        const size_t child = find_child_by_first_char(current, pattern[pos]);
        if (child == npos)
          return {};
 
        const size_t edge_start = nodes_[child].start;
        const size_t edge_end = nodes_[child].end;
 
        size_t k = 0;
        while (pos + k < pattern.size() and edge_start + k < edge_end
               and text_[edge_start + k] != text_.back()  // Never match sentinel
               and pattern[pos + k] == text_[edge_start + k])
          ++k;
 
        if (pos + k == pattern.size())
          {
            Array<size_t> matches;
            collect_leaf_suffixes(child, matches);
            suffix_structures_detail::sort_matches(matches);
            return matches;
          }
 
        if (edge_start + k == edge_end)
          {
            current = child;
            pos += k;
            continue;
          }
 
        return {};
      }
 
    return {};
  }
 
  [[nodiscard]] size_t node_count() const noexcept
  {
    return nodes_.size();
  }
 
  [[nodiscard]] size_t text_size() const noexcept
  {
    return original_size_;
  }
 
  [[nodiscard]] const Array<Node> &nodes() const noexcept
  {
    return nodes_;
  }
};
 
class Suffix_Automaton
{
public:
  struct State
  {
    std::array<int, 256> next; 
    int link = -1;             
    size_t len = 0;            
    size_t first_pos = 0;      
    bool terminal = false;     
    State()
    {
      next.fill(-1);
    }
  };
 
private:
  Array<State> states_;
  int last_ = 0;
 
  void mark_terminals()
  {
    int p = last_;
    while (p != -1)
      {
        states_[static_cast<size_t>(p)].terminal = true;
        p = states_[static_cast<size_t>(p)].link;
      }
  }
 
public:
  Suffix_Automaton()
  {
    clear();
  }
 
  void clear()
  {
    states_.empty();
    states_.append(State{});
    states_[0].link = -1;
    states_[0].len = 0;
    states_[0].first_pos = 0;
    states_[0].terminal = false;
    last_ = 0;
  }
 
  void extend(const char ch)
  {
    const auto c = static_cast<unsigned char>(ch);
 
    const int cur = static_cast<int>(states_.size());
    states_.append(State{});
    states_[static_cast<size_t>(cur)].len = states_[static_cast<size_t>(last_)].len + 1;
    states_[static_cast<size_t>(cur)].first_pos = states_[static_cast<size_t>(cur)].len - 1;
 
    int p = last_;
    while (p != -1 and states_[static_cast<size_t>(p)].next[c] == -1)
      {
        states_[static_cast<size_t>(p)].next[c] = cur;
        p = states_[static_cast<size_t>(p)].link;
      }
 
    if (p == -1)
      states_[static_cast<size_t>(cur)].link = 0;
    else
      {
        if (const int q = states_[static_cast<size_t>(p)].next[c];
            states_[static_cast<size_t>(p)].len + 1 == states_[static_cast<size_t>(q)].len)
          states_[static_cast<size_t>(cur)].link = q;
        else
          {
            const int clone = static_cast<int>(states_.size());
            // Copy q's state before append: append may reallocate states_,
            // invalidating any reference into the old buffer.
            const State clone_state = states_[static_cast<size_t>(q)];
            states_.append(clone_state);
            states_[static_cast<size_t>(clone)].len = states_[static_cast<size_t>(p)].len + 1;
            states_[static_cast<size_t>(clone)].terminal = false;
            // Enforce invariant: only root (state 0) may have link == -1
            if (states_[static_cast<size_t>(clone)].link == -1)
              states_[static_cast<size_t>(clone)].link = 0;
 
            while (p != -1 and states_[static_cast<size_t>(p)].next[c] == q)
              {
                states_[static_cast<size_t>(p)].next[c] = clone;
                p = states_[static_cast<size_t>(p)].link;
              }
 
            states_[static_cast<size_t>(q)].link = clone;
            states_[static_cast<size_t>(cur)].link = clone;
          }
      }
 
    last_ = cur;
  }
 
  void build(const std::string_view text)
  {
    clear();
    for (const char c : text)
      extend(c);
 
    mark_terminals();
  }
 
  [[nodiscard]] bool contains(const std::string_view pattern) const
  {
    int state = 0;
    for (const unsigned char c : pattern)
      {
        state = states_[static_cast<size_t>(state)].next[c];
        if (state == -1)
          return false;
      }
 
    return true;
  }
 
  [[nodiscard]] size_t distinct_substring_count() const
  {
    size_t total = 0;
    for (size_t i = 1; i < states_.size(); ++i)
      {
        const State &s = states_[i];
        const State &parent = states_[static_cast<size_t>(s.link)];
        total += s.len - parent.len;
      }
 
    return total;
  }
 
  [[nodiscard]] std::string longest_common_substring(const std::string_view other) const
  {
    int state = 0;
    size_t len = 0;
 
    size_t best_len = 0;
    size_t best_pos = 0;
 
    for (size_t i = 0; i < other.size(); ++i)
      {
        const auto c = static_cast<unsigned char>(other[i]);
 
        while (state != 0 and states_[static_cast<size_t>(state)].next[c] == -1)
          {
            state = states_[static_cast<size_t>(state)].link;
            len = states_[static_cast<size_t>(state)].len;
          }
 
        if (states_[static_cast<size_t>(state)].next[c] != -1)
          {
            state = states_[static_cast<size_t>(state)].next[c];
            ++len;
          }
        else
          {
            state = 0;
            len = 0;
          }
 
        if (len > best_len)
          {
            best_len = len;
            best_pos = i;
          }
      }
 
    if (best_len == 0)
      return "";
 
    return std::string(other.substr(best_pos + 1 - best_len, best_len));
  }
 
  [[nodiscard]] size_t state_count() const noexcept
  {
    return states_.size();
  }
 
  [[nodiscard]] const Array<State> &states() const noexcept
  {
    return states_;
  }
};
 
inline std::string longest_common_substring_sam(const std::string_view a, const std::string_view b)
{
  Suffix_Automaton sam;
  sam.build(a);
  return sam.longest_common_substring(b);
}
}  // namespace Aleph
 
#endif  // SUFFIX_STRUCTURES_H