dfs_pie_accumulate Subroutine

   recursive subroutine dfs_pie_accumulate(primary_atoms, primary_n_atoms, n_primaries, max_atoms, &
                                           clique, clique_size, current_atoms, n_current_atoms, &
                                           candidates, n_candidates, max_k_level, &
                                           atom_sets, coefficients, n_terms, max_terms, error)
      !! DFS helper: accumulate PIE coefficients for intersections
      integer, intent(in) :: primary_atoms(:, :)    !! Precomputed atom lists
      integer, intent(in) :: primary_n_atoms(:)     !! Atom counts
      integer, intent(in) :: n_primaries, max_atoms
      integer, intent(in) :: clique(:)              !! Current clique
      integer, intent(in) :: clique_size            !! Size of current clique
      integer, intent(in) :: current_atoms(:)       !! Atoms in current intersection
      integer, intent(in) :: n_current_atoms        !! Number of atoms in intersection
      integer, intent(in) :: candidates(:)          !! Candidate primaries
      integer, intent(in) :: n_candidates
      integer, intent(in) :: max_k_level
      integer, allocatable, intent(inout) :: atom_sets(:, :)
      integer, allocatable, intent(inout) :: coefficients(:)
      integer(int64), intent(inout) :: n_terms
      integer(int64), intent(inout) :: max_terms
      type(error_t), intent(inout) :: error         !! Error status

      integer :: sign, i, candidate_idx, candidate_pos
      integer(int64) :: term_idx
      integer, allocatable :: new_atoms(:), new_candidates(:)
      integer :: n_new_atoms, n_new_candidates
      logical :: found

      ! Skip empty intersections
      if (n_current_atoms == 0) return

      ! Compute PIE sign: (+1) for odd clique size, (-1) for even
      sign = merge(1, -1, mod(clique_size, 2) == 1)

      ! Find or create entry for this atom set
      found = .false.
      do term_idx = 1_int64, n_terms
         if (atom_sets_equal(atom_sets(:, term_idx), current_atoms, n_current_atoms)) then
            coefficients(term_idx) = coefficients(term_idx) + sign
            found = .true.
            exit
         end if
      end do

      if (.not. found) then
         ! New atom set
         if (n_terms >= max_terms) then
            call grow_pie_storage(atom_sets, coefficients, max_terms, max_atoms, error)
            if (error%has_error()) return
         end if
         n_terms = n_terms + 1_int64
         atom_sets(1:n_current_atoms, n_terms) = current_atoms(1:n_current_atoms)
         atom_sets(n_current_atoms + 1:, n_terms) = -1
         coefficients(n_terms) = sign
      end if

      ! Stop if we've reached maximum clique size
      if (clique_size >= max_k_level) return
      if (n_candidates == 0) return

      ! Try adding each candidate to the clique
      allocate (new_atoms(max_atoms))
      allocate (new_candidates(n_primaries))

      do i = 1, n_candidates
         candidate_idx = candidates(i)

         ! Compute intersection with this candidate
         call intersect_atom_lists(current_atoms, n_current_atoms, &
                                   primary_atoms(:, candidate_idx), primary_n_atoms(candidate_idx), &
                                   new_atoms, n_new_atoms)

         ! Skip if no intersection
         if (n_new_atoms == 0) cycle

         ! New candidates: must come after this one and overlap with new_atoms
         n_new_candidates = 0
         do candidate_pos = i + 1, n_candidates
            block
               integer :: test_candidate, test_n_intersect
               integer, allocatable :: test_intersect(:)
               test_candidate = candidates(candidate_pos)

               allocate (test_intersect(max_atoms))
               call intersect_atom_lists(new_atoms, n_new_atoms, &
                                         primary_atoms(:, test_candidate), primary_n_atoms(test_candidate), &
                                         test_intersect, test_n_intersect)

               if (test_n_intersect > 0) then
                  n_new_candidates = n_new_candidates + 1
                  new_candidates(n_new_candidates) = test_candidate
               end if
               deallocate (test_intersect)
            end block
         end do

         ! Recurse
         block
            integer :: new_clique(clique_size + 1)
            new_clique(1:clique_size) = clique(1:clique_size)
            new_clique(clique_size + 1) = candidate_idx

            call dfs_pie_accumulate(primary_atoms, primary_n_atoms, n_primaries, max_atoms, &
                                    new_clique, clique_size + 1, new_atoms, n_new_atoms, &
                                    new_candidates, n_new_candidates, max_k_level, &
                                    atom_sets, coefficients, n_terms, max_terms, error)
            if (error%has_error()) then
               deallocate (new_atoms, new_candidates)
               return
            end if
         end block
      end do

      deallocate (new_atoms, new_candidates)

   end subroutine dfs_pie_accumulate