#include <iostream>
#include <vector>
#include <array>
#include <tuple>
#include <algorithm>
#include <cmath>
#include <thread>
#include <numeric>
#include <memory>

// Requested constants
const int m = 4, n = 4;
const int cnt_monsters = 2;
const int antiverbosity = 100;

const int board_size = m * n;
const int cnt_states = std::pow(board_size, cnt_monsters) + .5;
int directions[5][2] = {{-1, 0}, {1, 0}, {0, -1}, {0, 1}, {0, 0}};
const int max_threads = 24;

const double monster_move_prob = std::pow(0.2, cnt_monsters);
const int threads_to_use = std::min(max_threads, cnt_states);

int get_state_index(const std::array<int, cnt_monsters>& state) {
  int index = 0;
  for (int i = 0; i < cnt_monsters; ++i) {
    index = index * board_size + state[i];
  }
  return index;
}

std::array<int, cnt_monsters> get_state_from_index(int index) {
  std::array<int, cnt_monsters> state;
  for (int i = cnt_monsters - 1; i >= 0; --i) {
    state[i] = index % board_size;
    index /= board_size;
  }
  return state;
}

void update_dp_by_state(int state_start, int state_end, const std::vector<std::vector<double>>& dp, std::vector<std::vector<double>>& new_dp, std::vector<std::vector<int>>& next_move) {
  for (int state_idx = state_start; state_idx < state_end; ++state_idx) {
    // OPTIMIZATION: Calculate monster outcomes once, as they only depend on state_idx.
    std::vector<std::array<int, cnt_monsters>> monster_outcomes;
    monster_outcomes.push_back(get_state_from_index(state_idx));

    for (int i = 0; i < cnt_monsters; ++i) {
      std::vector<std::array<int, cnt_monsters>> generated_states;
      generated_states.reserve(monster_outcomes.size() * 5);
      for (const auto& current_state : monster_outcomes) {
        int monster_x = current_state[i] / n, monster_y = current_state[i] % n;
        for (int monster_d = 0; monster_d < 5; ++monster_d) {
          auto [mdx, mdy] = directions[monster_d];
          auto next_state = current_state;
          int new_monster_x = monster_x + mdx;
          int new_monster_y = monster_y + mdy;

          if (0 <= new_monster_x && new_monster_x < m && 0 <= new_monster_y && new_monster_y < n) {
            next_state[i] = new_monster_x * n + new_monster_y;
          }
          generated_states.push_back(next_state);
        }
      }
      monster_outcomes.swap(generated_states);
    }

    // Inner loop over man positions.
    for (int man_pos = 0; man_pos < board_size; ++man_pos) {
      if (dp[man_pos][state_idx] == -1) continue;

      int mx = man_pos / n, my = man_pos % n;
      std::array<double, 5> expected_scores{};
      for (int d = 0; d < 5; ++d) {
        auto [dx, dy] = directions[d];
        int new_mx = mx + dx, new_my = my + dy;

        if (new_mx < 0 || new_mx >= m || new_my < 0 || new_my >= n) {
          expected_scores[d] = -1.0;
          continue;
        }
        int new_man_pos = new_mx * n + new_my;

        double total_score = 0;
        for (const auto& new_monster_state : monster_outcomes) {
          int new_state_idx = get_state_index(new_monster_state);
          if (dp[new_man_pos][new_state_idx] != -1) {
            total_score += monster_move_prob * (1 + dp[new_man_pos][new_state_idx]);
          }
        }
        expected_scores[d] = total_score;
      }
      int best_move = std::distance(expected_scores.begin(), std::max_element(expected_scores.begin(), expected_scores.end()));
      new_dp[man_pos][state_idx] = expected_scores[best_move];
      next_move[man_pos][state_idx] = best_move;
    }
  }
}

void compute_divergence_by_state(int state_start, int state_end,
  const std::vector<std::vector<double>>& old_dp, const std::vector<std::vector<double>>& current_dp,
  double& partial_divergence, double& partial_avg_sum, int& partial_valid_count) {

  double div = 0, avg_sum = 0;
  int count = 0;
  for (int state_idx = state_start; state_idx < state_end; ++state_idx) {
    for (int man_pos = 0; man_pos < board_size; ++man_pos) {
      if (old_dp[man_pos][state_idx] != -1) {
        double val1 = old_dp[man_pos][state_idx];
        double val2 = current_dp[man_pos][state_idx];
        div += (val1 > val2) ? (val1 - val2) : (val2 - val1);
        avg_sum += current_dp[man_pos][state_idx];
        count++;
      }
    }
  }
  partial_divergence = div;
  partial_avg_sum = avg_sum;
  partial_valid_count = count;
}

int main() {
  std::cout << "Board: " << m << "x" << n << ", Monsters: " << cnt_monsters << ", Threads: " << threads_to_use << std::endl;
  std::cout << "State space size: " << board_size * cnt_states << std::endl;

  auto dp = std::make_unique<std::vector<std::vector<double>>>(board_size, std::vector<double>(cnt_states, 0));
  auto new_dp = std::make_unique<std::vector<std::vector<double>>>(board_size, std::vector<double>(cnt_states, 0));
  auto next_move = std::make_unique<std::vector<std::vector<int>>>(board_size, std::vector<int>(cnt_states, 0));

  for (int man_pos = 0; man_pos < board_size; ++man_pos) {
    for (int state_idx = 0; state_idx < cnt_states; ++state_idx) {
      auto monster_state = get_state_from_index(state_idx);
      for (int pos : monster_state) {
        if (pos == man_pos) {
          (*dp)[man_pos][state_idx] = -1;
          (*new_dp)[man_pos][state_idx] = -1;
          break;
        }
      }
    }
  }

  for (int t = 0;; ++t) {
    if (t % antiverbosity == 0) {
      std::cout << "Iteration " << t << std::endl;
    }

    int chunk_size = cnt_states / threads_to_use;
    std::vector<std::thread> threads;
    for (int i = 0; i < threads_to_use; ++i) {
      int start = i * chunk_size;
      int end = (i == threads_to_use - 1) ? cnt_states : (i + 1) * chunk_size;
      threads.emplace_back(update_dp_by_state, start, end, std::cref(*dp), std::ref(*new_dp), std::ref(*next_move));
    }
    for (auto& thread : threads) {
      thread.join();
    }

    std::vector<std::thread> div_threads;
    std::vector<double> partial_divergences(threads_to_use);
    std::vector<double> partial_avg_sums(threads_to_use);
    std::vector<int> partial_valid_counts(threads_to_use);

    for (int i = 0; i < threads_to_use; ++i) {
      int start = i * chunk_size;
      int end = (i == threads_to_use - 1) ? cnt_states : (i + 1) * chunk_size;
      div_threads.emplace_back(compute_divergence_by_state, start, end,
        std::cref(*dp), std::cref(*new_dp),
        std::ref(partial_divergences[i]),
        std::ref(partial_avg_sums[i]),
        std::ref(partial_valid_counts[i]));
    }
    for (auto& thread : div_threads) {
      thread.join();
    }

    double divergence = std::accumulate(partial_divergences.begin(), partial_divergences.end(), 0.0);
    double avg_sum = std::accumulate(partial_avg_sums.begin(), partial_avg_sums.end(), 0.0);
    int valid_states_count = std::accumulate(partial_valid_counts.begin(), partial_valid_counts.end(), 0);
    double avg = (valid_states_count > 0) ? avg_sum / valid_states_count : 0;

    if (t % antiverbosity == 0) {
      std::cout << "Divergence " << divergence << ", average survival time " << avg << std::endl;
    }
    if (divergence < 1e-1) break;

    dp.swap(new_dp);
  }

  // Example: man at (3,3) [index 15], monsters at (0,0) [0] and (0,1) [4]
  // state_index = 0 * 16 + 4 = 4
  double value = (*dp)[15][4];
  std::cout << "Expected survival time from man at (3,3), monsters at (0,0) & (0,1): " << value << std::endl;

  return 0;
}


rewrite to PURE C, and also, remove excessive comments and one-line functions