.. _program_listing_file_src_cpp_parallel_factorizer.cpp:

Program Listing for File parallel_factorizer.cpp
================================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_cpp_parallel_factorizer.cpp>` (``src/cpp/parallel_factorizer.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "parallel_factorizer.hpp"
   #include "factorizer_helpers.hpp"
   #include <algorithm>
   #include <iostream>
   #include <string>
   #include <chrono>
   
   namespace fs = std::filesystem;
   namespace noLZSS {
   // Helper functions lcp() and next_leaf() are now in factorizer_helpers.hpp
   
   std::string ParallelFactorizer::create_temp_file_path(size_t thread_id) {
       auto timestamp = std::chrono::system_clock::now().time_since_epoch().count();
       std::string temp_dir = fs::temp_directory_path().string();
       return temp_dir + "/noLZSS_temp_" + std::to_string(timestamp) + "_" + 
              std::to_string(thread_id) + ".bin";
   }
   
   size_t ParallelFactorizer::parallel_factorize(std::string_view text, const std::string& output_path, 
                                              size_t num_threads, size_t start_pos) {
       if (text.empty()) return 0;
       
       // Ensure start_pos is within bounds
       if (start_pos >= text.length()) {
           throw std::invalid_argument("start_pos must be less than text length");
       }
       
       // Determine optimal thread count if not specified
       if (num_threads == 0) {
           // Calculate maximum useful threads based on remaining text size
           size_t remaining_length = text.length() - start_pos;
           size_t max_useful_threads = remaining_length / MIN_CHARS_PER_THREAD;
           
           // Use available hardware threads, but don't exceed what's useful for this input size
           size_t hardware_threads = std::thread::hardware_concurrency();
           num_threads = std::min(hardware_threads, max_useful_threads);
           
           // Ensure at least one thread
           num_threads = std::max(1UL, num_threads);
       }
       
       // Create suffix tree for all threads to use
       std::string tmp(text);
       cst_t cst;
       construct_im(cst, tmp, 1);
       
       // Create RMQ support once for all threads
       sdsl::rmq_succinct_sct<> rmq(&cst.csa);
       
       // Create contexts for each thread
       std::vector<ThreadContext> contexts(num_threads);
       
       // Divide work among threads, starting from start_pos
       const size_t remaining_length = text.length() - start_pos;
       const size_t chunk_size = remaining_length / num_threads;
       for (size_t i = 0; i < num_threads; ++i) {
           contexts[i].thread_id = i;
           contexts[i].start_pos = start_pos + (i * chunk_size);
           contexts[i].end_pos = (i + 1 < num_threads) ? (start_pos + (i + 1) * chunk_size) : text.length();
           contexts[i].text_length = text.length();
           // Thread 0 writes directly to output file; others use temp files
           contexts[i].temp_file_path = (i == 0) ? output_path : create_temp_file_path(i);
           contexts[i].is_last_thread = (i == num_threads - 1);
       }
       
       // Create mutexes for file access
       std::vector<std::mutex> file_mutexes(num_threads);
       
       // Track exceptions from threads
       std::vector<std::exception_ptr> thread_exceptions(num_threads, nullptr);
       std::mutex exception_mutex;
       
       // Create and start worker threads with exception handling
       std::vector<std::thread> threads;
       for (size_t i = 0; i < num_threads; ++i) {
           threads.emplace_back([this, &cst, &rmq, &contexts, &file_mutexes, &thread_exceptions, &exception_mutex, i]() {
               try {
                   factorize_thread(cst, rmq, contexts[i], contexts, file_mutexes);
               } catch (...) {
                   // Capture exception for the main thread to handle
                   std::lock_guard<std::mutex> lock(exception_mutex);
                   thread_exceptions[i] = std::current_exception();
               }
           });
       }
       
       // Wait for all threads to complete
       for (auto& t : threads) {
           if (t.joinable()) t.join();
       }
       
       // Check for exceptions from threads
       for (size_t i = 0; i < num_threads; ++i) {
           if (thread_exceptions[i]) {
               // Cleanup temporary files before rethrowing
               cleanup_temp_files(contexts);
               std::rethrow_exception(thread_exceptions[i]);
           }
       }
       
       // Merge temporary files and create final output
       size_t total_factors = merge_temp_files(output_path, contexts);
       
       // Cleanup temporary files
       cleanup_temp_files(contexts);
       
       return total_factors;
   }
   
   void ParallelFactorizer::factorize_thread(const cst_t& cst, const sdsl::rmq_succinct_sct<>& rmq,
                                           ThreadContext& ctx,
                                           std::vector<ThreadContext>& all_contexts,
                                           std::vector<std::mutex>& file_mutexes) {
       // Initialize factorization at our starting position
       auto lambda = cst.select_leaf(cst.csa.isa[ctx.start_pos] + 1);
       size_t lambda_node_depth = cst.node_depth(lambda);
       size_t lambda_sufnum = ctx.start_pos;
       
       // Create or truncate output file
       {
           std::lock_guard<std::mutex> lock(file_mutexes[ctx.thread_id]);
           std::ofstream ofs(ctx.temp_file_path, std::ios::binary | std::ios::trunc);
       }
       
       // Track the next thread for convergence checking
       ThreadContext* next_ctx = nullptr;
       if (!ctx.is_last_thread && ctx.thread_id + 1 < all_contexts.size()) {
           next_ctx = &all_contexts[ctx.thread_id + 1];
       }
       
       // Main factorization loop
       while (lambda_sufnum < ctx.text_length) {
           // Compute current factor
           size_t d = 1;
           size_t l = 1;
           cst_t::node_type v;
           size_t v_min_leaf_sufnum = 0;
           size_t u_min_leaf_sufnum = 0;
           Factor current_factor;
           
           // Factor computation logic (similar to original nolzss)
           while (true) {
               v = cst.bp_support.level_anc(lambda, lambda_node_depth - d);
               v_min_leaf_sufnum = cst.csa[rmq(cst.lb(v), cst.rb(v))];
               l = cst.depth(v);
               
               if (v_min_leaf_sufnum + l - 1 < lambda_sufnum) {
                   u_min_leaf_sufnum = v_min_leaf_sufnum;
                   ++d;
                   continue;
               }
               
               auto u = cst.parent(v);
               auto u_depth = cst.depth(u);
               
               if (v_min_leaf_sufnum == lambda_sufnum) {
                   if (u == cst.root()) {
                       l = 1;
                       current_factor = Factor{static_cast<uint64_t>(lambda_sufnum), 
                                            static_cast<uint64_t>(l), 
                                            static_cast<uint64_t>(lambda_sufnum)};
                   } else {
                       l = u_depth;
                       current_factor = Factor{static_cast<uint64_t>(lambda_sufnum), 
                                            static_cast<uint64_t>(l), 
                                            static_cast<uint64_t>(u_min_leaf_sufnum)};
                   }
               } else {
                   l = std::min(lcp(cst, lambda_sufnum, v_min_leaf_sufnum),
                              (lambda_sufnum - v_min_leaf_sufnum));
                   if (l <= u_depth) {
                       l = u_depth;
                       current_factor = Factor{static_cast<uint64_t>(lambda_sufnum), 
                                            static_cast<uint64_t>(l), 
                                            static_cast<uint64_t>(u_min_leaf_sufnum)};
                   } else {
                       current_factor = Factor{static_cast<uint64_t>(lambda_sufnum), 
                                            static_cast<uint64_t>(l), 
                                            static_cast<uint64_t>(v_min_leaf_sufnum)};
                   }
               }
               break;
           }
           
           // Write the factor to temporary file
           write_factor(current_factor, ctx.temp_file_path, file_mutexes[ctx.thread_id]);
           
           // Check for convergence when factor extends into next thread's region
           // Use lambda_sufnum + l (exclusive end) to check against next thread's start
           if (next_ctx && lambda_sufnum + l >= next_ctx->start_pos) {
               size_t current_end = lambda_sufnum + l;  // Exclusive end position
               if (check_convergence(current_end, *next_ctx, file_mutexes[next_ctx->thread_id])) {
                   break;  // Convergence detected
               }
           }
           
           // Advance to next position
           lambda = next_leaf(cst, lambda, l);
           lambda_node_depth = cst.node_depth(lambda);
           lambda_sufnum = cst.sn(lambda);
       }
   }
   
   void ParallelFactorizer::factorize_dna_w_rc_thread(const cst_t& cst,
                                                     const sdsl::rmq_succinct_sct<>& rmqF,
                                                     const sdsl::rmq_succinct_sct<>& rmqRcEnd,
                                                     const sdsl::int_vector<64>& fwd_starts,
                                                     const sdsl::int_vector<64>& rc_ends,
                                                     uint64_t INF,
                                                     size_t N,
                                                     ThreadContext& ctx,
                                                     std::vector<ThreadContext>& all_contexts,
                                                     std::vector<std::mutex>& file_mutexes) {
       // Initialize to the leaf of suffix starting at position ctx.start_pos
       auto lambda = cst.select_leaf(cst.csa.isa[ctx.start_pos] + 1);
       size_t lambda_node_depth = cst.node_depth(lambda);
       size_t i = cst.sn(lambda); // Current position in text
       
       // Create or truncate output file
       {
           std::lock_guard<std::mutex> lock(file_mutexes[ctx.thread_id]);
           std::ofstream ofs(ctx.temp_file_path, std::ios::binary | std::ios::trunc);
       }
       
       // Track the next thread for convergence checking
       ThreadContext* next_ctx = nullptr;
       if (!ctx.is_last_thread && ctx.thread_id + 1 < all_contexts.size()) {
           next_ctx = &all_contexts[ctx.thread_id + 1];
       }
       
       // Main factorization loop - only process the original sequence (0 to N)
       while (i < N) {
           // At factor start i (0-based in T), walk up ancestors and track best FWD and RC candidates independently
           
           // Track best FWD candidate
           bool   have_fwd       = false;
           size_t best_fwd_start = 0;  // start in T (for FWD)
           size_t best_fwd_depth = 0;  // tree depth of best FWD candidate
           
           // Track best RC candidate
           bool   have_rc        = false;
           size_t best_rc_end    = 0;  // end in T (for RC)
           size_t best_rc_posS   = 0;  // pos in S where RC candidate suffix starts (for LCP)
           size_t best_rc_depth  = 0;  // tree depth of best RC candidate
   
           // Walk from leaf to root via level_anc
           for (size_t step = 1; step <= lambda_node_depth; ++step) {
               auto v = cst.bp_support.level_anc(lambda, lambda_node_depth - step);
               size_t ell = cst.depth(v);
               if (ell == 0) break; // reached root
   
               auto lb = cst.lb(v), rb = cst.rb(v);
   
               // Forward candidate (min start in T within v's interval)
               size_t kF = rmqF(lb, rb);
               uint64_t jF = fwd_starts[kF];
               bool okF = (jF != INF) && (jF + ell - 1 < i); // non-overlap: endF <= i-1
   
               // RC candidate (min END in T within v's interval; monotone with depth)
               size_t kR = rmqRcEnd(lb, rb);
               uint64_t endRC = rc_ends[kR];
               bool okR = (endRC != INF) && (endRC < i); // endRC <= i-1
   
               if (!okF && !okR) {
                   // deeper nodes can only increase jF and the minimal RC end
                   // -> non-overlap won't become true again for either; stop
                   break;
               }
   
               // Update best FWD candidate independently
               if (okF) {
                   if (ell > best_fwd_depth ||
                       (ell == best_fwd_depth && (jF + ell - 1) < (best_fwd_start + best_fwd_depth - 1))) {
                       have_fwd = true;
                       best_fwd_start = jF;
                       best_fwd_depth = ell;
                   }
               }
               
               // Update best RC candidate independently
               if (okR) {
                   size_t posS_R = cst.csa[kR]; // suffix position in S for LCP
                   if (ell > best_rc_depth ||
                       (ell == best_rc_depth && (endRC < best_rc_end))) {
                       have_rc = true;
                       best_rc_end = endRC;
                       best_rc_posS = posS_R;
                       best_rc_depth = ell;
                   }
               }
           }
   
           size_t emit_len = 1;
           uint64_t emit_ref = i; // default for literal
           
           if (!have_fwd && !have_rc) {
               // No previous occurrence (FWD nor RC) — literal of length 1
               Factor f{static_cast<uint64_t>(i), static_cast<uint64_t>(emit_len), static_cast<uint64_t>(emit_ref)};
               write_factor(f, ctx.temp_file_path, file_mutexes[ctx.thread_id]);
           } else {
               // Compute true lengths for both candidates
               size_t fwd_true_len = 0;
               size_t rc_true_len = 0;
               
               if (have_fwd) {
                   size_t cap = i - best_fwd_start; // non-overlap cap
                   size_t L = lcp(cst, i, best_fwd_start);
                   fwd_true_len = std::min(L, cap);
               }
               
               if (have_rc) {
                   rc_true_len = lcp(cst, i, best_rc_posS);
               }
               
               // Choose based on TRUE length, with FWD preference at equal length.
               // A literal (self-reference, length 1) is always available as a fallback.
               // Only use RC if it provides strictly longer match than both FWD and literal.
               bool use_fwd = false;
               bool use_literal = false;
               
               if (have_fwd && fwd_true_len >= 1) {
                   // Have a real forward match
                   if (have_rc && rc_true_len > fwd_true_len) {
                       use_fwd = false;  // RC is strictly longer
                   } else {
                       use_fwd = true;   // FWD wins (including ties)
                   }
               } else {
                   // No real forward match; compare RC against literal (length 1)
                   if (have_rc && rc_true_len > 1) {
                       use_fwd = false;  // RC is strictly longer than literal
                   } else {
                       use_literal = true;  // Use literal (self-reference)
                   }
               }
               
               if (use_literal) {
                   // Emit literal
                   emit_len = 1;
                   emit_ref = static_cast<uint64_t>(i);
               } else if (use_fwd) {
                   emit_len = fwd_true_len;
                   emit_ref = static_cast<uint64_t>(best_fwd_start);
               } else {
                   emit_len = rc_true_len;
                   size_t start_pos_val = best_rc_end - emit_len + 1;  // start = end - len + 1
                   emit_ref = RC_MASK | static_cast<uint64_t>(start_pos_val); // start-anchored + RC flag
               }
               
               // Safety: ensure progress
               if (emit_len == 0) {
                   throw std::runtime_error("emit_len must be positive to ensure factorization progress");
               }
   
               Factor f{static_cast<uint64_t>(i), static_cast<uint64_t>(emit_len), emit_ref};
               write_factor(f, ctx.temp_file_path, file_mutexes[ctx.thread_id]);
           }
           
           // Check for convergence
           if (next_ctx && i + emit_len >= next_ctx->start_pos) {
               size_t current_end = i + emit_len;
               if (check_convergence(current_end, *next_ctx, file_mutexes[next_ctx->thread_id])) {
                   break; // Convergence detected
               }
           }
           
           // Advance to next position
           lambda = next_leaf(cst, lambda, emit_len);
           lambda_node_depth = cst.node_depth(lambda);
           i = cst.sn(lambda);
       }
   }
   
   std::optional<Factor> ParallelFactorizer::read_factor_at_index(const std::string& file_path, 
                                                                   size_t factor_index) {
       std::ifstream ifs(file_path, std::ios::binary);
       if (!ifs) {
           return std::nullopt;
       }
       
       // Seek to the position of the requested factor
       ifs.seekg(factor_index * sizeof(Factor), std::ios::beg);
       if (!ifs) {
           return std::nullopt; // Index out of bounds
       }
       
       // Read the factor
       Factor factor;
       ifs.read(reinterpret_cast<char*>(&factor), sizeof(Factor));
       if (!ifs || ifs.gcount() != sizeof(Factor)) {
           return std::nullopt; // Failed to read complete factor
       }
       
       return factor;
   }
   
   bool ParallelFactorizer::check_convergence(size_t current_end, ThreadContext& next_ctx,
                                           std::mutex& next_file_mutex) {
       // If we have a cached factor, check it first
       if (next_ctx.last_read_factor.has_value()) {
           const auto& cached_factor = next_ctx.last_read_factor.value();
           size_t cached_end = cached_factor.start + cached_factor.length;
           
           if (cached_end > current_end) {
               // The cached factor ends beyond our current position,
               // so we haven't reached convergence yet
               return false;
           }
           
           if (cached_factor.start == current_end) {
               // Found convergence: next thread's factor starts exactly where we end
               return true;
           }
           
           // cached_end <= current_end but cached_factor.start < current_end
           // Need to read the next factor
       }
       
       // Read factors one at a time until we find convergence or go past current_end
       while (true) {
           std::lock_guard<std::mutex> lock(next_file_mutex);
           
           auto factor_opt = read_factor_at_index(next_ctx.temp_file_path, 
                                                  next_ctx.next_thread_factor_index);
           
           if (!factor_opt.has_value()) {
               // No more factors in next thread's file - convergence not found yet
               return false;
           }
           
           Factor factor = factor_opt.value();
           next_ctx.last_read_factor = factor;
           next_ctx.next_thread_factor_index++;
           
           if (factor.start == current_end) {
               // Found convergence
               return true;
           }
           
           if (factor.start > current_end) {
               // Next thread's factor starts beyond our current position
               // We haven't reached convergence yet
               return false;
           }
           
           // factor.start < current_end, keep reading next factors
       }
   }
   
   void ParallelFactorizer::write_factor(const Factor& factor, const std::string& file_path, 
                                      std::mutex& file_mutex) {
       std::lock_guard<std::mutex> lock(file_mutex);
       
       std::ofstream ofs(file_path, std::ios::binary | std::ios::app);
       if (!ofs) {
           throw std::runtime_error("Cannot open temporary file for writing: " + file_path);
       }
       
       ofs.write(reinterpret_cast<const char*>(&factor), sizeof(Factor));
   }
   
   std::optional<Factor> ParallelFactorizer::read_factor_at(const std::string& file_path, 
                                                        size_t index, 
                                                        std::mutex& file_mutex) {
       std::lock_guard<std::mutex> lock(file_mutex);
       
       std::ifstream ifs(file_path, std::ios::binary);
       if (!ifs) {
           return std::nullopt;
       }
       
       // Seek to the position of the requested factor
       ifs.seekg(index * sizeof(Factor));
       
       Factor factor;
       if (ifs.read(reinterpret_cast<char*>(&factor), sizeof(Factor))) {
           return factor;
       }
       
       return std::nullopt;
   }
   
   size_t ParallelFactorizer::merge_temp_files(const std::string& output_path,
                                            std::vector<ThreadContext>& contexts) {
       // Thread 0 already wrote to output_path, so we open in append mode
       // to add factors from other threads
       
       size_t total_factors = 0;
       uint64_t total_length = 0;
       size_t current_position = 0;  // Track the current end position in the merged output
       size_t text_length = contexts.empty() ? 0 : contexts[0].text_length;
       std::optional<Factor> last_written_factor;
       
       // First, read thread 0's factors to find the last one and count them
       {
           std::ifstream ifs(output_path, std::ios::binary);
           if (!ifs) {
               throw std::runtime_error("Cannot open output file for reading: " + output_path);
           }
           
           Factor factor;
           while (ifs.read(reinterpret_cast<char*>(&factor), sizeof(Factor))) {
               total_factors++;
               total_length += factor.length;
               last_written_factor = factor;
               
               size_t factor_end = factor.start + factor.length;
               if (factor_end > current_position) {
                   current_position = factor_end;
               }
           }
       }
       
       // Process remaining threads if there are multiple threads and sequence isn't complete
       if (contexts.size() > 1 && current_position < text_length) {
           // Open output file in append mode for adding factors from other threads
           std::ofstream ofs(output_path, std::ios::binary | std::ios::app);
           if (!ofs) {
               throw std::runtime_error("Cannot open output file for appending: " + output_path);
           }
           
           // Process remaining threads (i >= 1)
           for (size_t i = 1; i < contexts.size(); i++) {
               // Check if we've already covered the entire sequence
               if (current_position >= text_length) {
                   break;  // Stop merging - sequence is complete
               }
               
               std::ifstream ifs(contexts[i].temp_file_path, std::ios::binary);
               if (!ifs) continue;
               
               Factor factor;
               bool found_convergence = false;
               
               // For subsequent threads, skip factors until we find convergence point
               // Convergence: last_written_factor.end == current_factor.start
               
               if (!last_written_factor.has_value()) {
                   // No last written factor - shouldn't happen, but handle gracefully
                   continue;
               }
               
               size_t last_end = last_written_factor->start + last_written_factor->length;
               
               // Read factors and look for convergence
               while (ifs.read(reinterpret_cast<char*>(&factor), sizeof(Factor))) {
                   if (!found_convergence) {
                       // Still looking for convergence point
                       if (factor.start == last_end) {
                           // Found convergence! Start copying from this factor onwards
                           found_convergence = true;
                       } else {
                           // Skip this factor - it's before convergence
                           continue;
                       }
                   }
                   
                   // Write factor to output (either we found convergence or we're continuing)
                   ofs.write(reinterpret_cast<const char*>(&factor), sizeof(Factor));
                   total_factors++;
                   total_length += factor.length;
                   last_written_factor = factor;
                   
                   size_t factor_end = factor.start + factor.length;
                   if (factor_end > current_position) {
                       current_position = factor_end;
                   }
                   last_end = factor_end;  // Update for next iteration
                   
                   // Stop if we've covered the entire sequence
                   if (current_position >= text_length) {
                       break;
                   }
               }
               
               // If we've covered the entire sequence, stop processing more threads
               if (current_position >= text_length) {
                   break;
               }
           }
       }
       
       // Write footer at the end (v2 format) - single location for all cases
       std::ofstream ofs(output_path, std::ios::binary | std::ios::app);
       if (!ofs) {
           throw std::runtime_error("Cannot open output file for appending footer: " + output_path);
       }
       
       FactorFileFooter footer;
       footer.num_factors = total_factors;
       footer.num_sequences = 0;  // Unknown for general (non-FASTA) factorization
       footer.num_sentinels = 0;  // No sentinels for general factorization
       footer.footer_size = sizeof(FactorFileFooter);
       footer.total_length = total_length;
       
       ofs.write(reinterpret_cast<const char*>(&footer), sizeof(footer));
       
       return total_factors;
   }
   
   void ParallelFactorizer::cleanup_temp_files(const std::vector<ThreadContext>& contexts) {
       for (const auto& ctx : contexts) {
           // Skip thread 0 - it wrote directly to the output file
           if (ctx.thread_id == 0) continue;
           
           try {
               fs::remove(ctx.temp_file_path);
           } catch (const std::exception& e) {
               std::cerr << "Warning: Failed to remove temporary file: " << ctx.temp_file_path
                        << " - " << e.what() << std::endl;
           }
       }
   }
   
   size_t ParallelFactorizer::parallel_factorize_file(const std::string& input_path, 
                                                   const std::string& output_path,
                                                   size_t num_threads,
                                                   size_t start_pos) {
       // Read the file content
       std::ifstream is(input_path, std::ios::binary);
       if (!is) {
           throw std::runtime_error("Cannot open input file: " + input_path);
       }
       
       std::string data((std::istreambuf_iterator<char>(is)), std::istreambuf_iterator<char>());
       return parallel_factorize(data, output_path, num_threads, start_pos);
   }
   
   size_t ParallelFactorizer::parallel_factorize_multiple_dna_w_rc(const std::string& prepared_string,
                                                                   size_t original_length,
                                                                   const std::string& output_path,
                                                                   size_t num_threads,
                                                                   size_t start_pos) {
       if (prepared_string.empty() || original_length == 0) return 0;
       
       const std::string& S = prepared_string;
       
       // Validate minimum string length to avoid underflow and SDSL errors
       // Format for w/ RC: T1!T2@...Tn$rt(Tn)%...rt(T1)&
       // Minimum valid input: "A$A&" (4 chars: sequence + sentinel + reverse_complement + sentinel)
       // - Position 0: Single nucleotide (e.g., 'A')
       // - Position 1: Sentinel separating original from RC (e.g., '$')
       // - Position 2: Reverse complement of the nucleotide (e.g., 'A' is revcomp of 'T')
       // - Position 3: Final sentinel (e.g., '&')
       if (S.size() < 4) {
           std::cerr << "Warning: Input string too short for factorization with reverse complement (size=" 
                     << S.size() << "). Returning 0 factors." << std::endl;
           return 0;
       }
   
       const size_t N = (S.size() / 2) - 1;
       if (N == 0) {
           std::cerr << "Warning: Computed N=0 from input size=" << S.size() 
                     << ". Returning 0 factors." << std::endl;
           return 0;
       }
       
       // Validate start_pos
       if (start_pos >= N) {
           throw std::invalid_argument("start_pos must be less than the original sequence length");
       }
       
       // For DNA w/ RC, we can only parallelize the original sequence part (start_pos to N)
       // The algorithm needs the full prepared string with RC for suffix tree
       
       // Determine optimal thread count if not specified
       if (num_threads == 0) {
           size_t remaining_length = N - start_pos;
           size_t max_useful_threads = remaining_length / MIN_CHARS_PER_THREAD;
           size_t hardware_threads = std::thread::hardware_concurrency();
           num_threads = std::min(hardware_threads, max_useful_threads);
           num_threads = std::max(1UL, num_threads);
       }
       
       // Build CST over the prepared string (T + sentinel + RC(T) + sentinel)
       cst_t cst;
       construct_im(cst, S, 1);
       
       // Build RMQ structures for DNA w/ RC algorithm
       const uint64_t INF = std::numeric_limits<uint64_t>::max() / 2ULL;
       sdsl::int_vector<64> fwd_starts(cst.csa.size(), INF);
       sdsl::int_vector<64> rc_ends(cst.csa.size(), INF);
       
       const size_t T_end = N;           // end of original (exclusive)
       const size_t R_beg = N + 1;       // first char of rc (skip sentinel at position N)
       const size_t R_end = S.size() - 1;   // exclude final sentinel
   
       for (size_t k = 0; k < cst.csa.size(); ++k) {
           size_t posS = cst.csa[k];
           if (posS < T_end) {
               // suffix starts in T
               fwd_starts[k] = posS;     // 0-based start in T
           } else if (posS >= R_beg && posS < R_end) {
               // suffix starts in R
               size_t jR0   = posS - R_beg;         // 0-based start in R
               size_t endT0 = N - jR0 - 1;          // mapped end in T (0-based)
               rc_ends[k] = endT0;
           }
       }
       
       sdsl::rmq_succinct_sct<> rmqF(&fwd_starts);
       sdsl::rmq_succinct_sct<> rmqRcEnd(&rc_ends);
       
       // Create contexts for each thread - divide only from start_pos to N
       std::vector<ThreadContext> contexts(num_threads);
       const size_t remaining_length = N - start_pos;
       const size_t chunk_size = remaining_length / num_threads;
       
       for (size_t i = 0; i < num_threads; ++i) {
           contexts[i].thread_id = i;
           contexts[i].start_pos = start_pos + (i * chunk_size);
           contexts[i].end_pos = (i + 1 < num_threads) ? (start_pos + (i + 1) * chunk_size) : N;
           contexts[i].text_length = N; // Only the original sequence length
           contexts[i].temp_file_path = (i == 0) ? output_path : create_temp_file_path(i);
           contexts[i].is_last_thread = (i == num_threads - 1);
       }
       
       // Create mutexes for file access
       std::vector<std::mutex> file_mutexes(num_threads);
       
       // Track exceptions from threads
       std::vector<std::exception_ptr> thread_exceptions(num_threads, nullptr);
       std::mutex exception_mutex;
       
       // Lambda to pass RMQ structures to threads with exception handling
       auto factorize_dna_thread = [&](ThreadContext& ctx) {
           try {
               factorize_dna_w_rc_thread(cst, rmqF, rmqRcEnd, fwd_starts, rc_ends, INF, N,
                                        ctx, contexts, file_mutexes);
           } catch (...) {
               // Capture exception for the main thread to handle
               std::lock_guard<std::mutex> lock(exception_mutex);
               thread_exceptions[ctx.thread_id] = std::current_exception();
           }
       };
       
       // Create and start worker threads
       std::vector<std::thread> threads;
       for (size_t i = 0; i < num_threads; ++i) {
           threads.emplace_back(factorize_dna_thread, std::ref(contexts[i]));
       }
       
       // Wait for all threads to complete
       for (auto& t : threads) {
           if (t.joinable()) t.join();
       }
       
       // Check for exceptions from threads
       for (size_t i = 0; i < num_threads; ++i) {
           if (thread_exceptions[i]) {
               // Cleanup temporary files before rethrowing
               cleanup_temp_files(contexts);
               std::rethrow_exception(thread_exceptions[i]);
           }
       }
       
       // Merge temporary files and create final output
       size_t total_factors = merge_temp_files(output_path, contexts);
       
       // Cleanup temporary files
       cleanup_temp_files(contexts);
       
       return total_factors;
   }
   
   size_t ParallelFactorizer::parallel_factorize_dna_w_rc(std::string_view text, 
                                                       const std::string& output_path,
                                                       size_t num_threads,
                                                       size_t start_pos) {
       if (text.empty()) return 0;
       
       // Prepare the DNA sequence with reverse complement
       std::string text_str(text);
       PreparedSequenceResult prep_result = prepare_multiple_dna_sequences_w_rc({text_str});
       
       // Call the core function with the prepared string
       return parallel_factorize_multiple_dna_w_rc(prep_result.prepared_string,
                                                  prep_result.original_length,
                                                  output_path,
                                                  num_threads,
                                                  start_pos);
   }
   
   size_t ParallelFactorizer::parallel_factorize_file_dna_w_rc(const std::string& input_path, 
                                                            const std::string& output_path,
                                                            size_t num_threads) {
       std::ifstream is(input_path, std::ios::binary);
       if (!is) {
           throw std::runtime_error("Cannot open input file: " + input_path);
       }
       
       std::string data((std::istreambuf_iterator<char>(is)), std::istreambuf_iterator<char>());
       return parallel_factorize_dna_w_rc(data, output_path, num_threads);
   }
   
   } // namespace noLZSS