Sudoku solver hangs for some games

I have written an implementation of Donald Knuth's Algorithm X for solving exact cover problems and applied it to Sudoku for the purpose of solving Project Euler Problem 96. My code is a translation from Python into C++ of Ali Assaf's implementation of the same algorithm.

My version is solves most of the grids contained in this text file, but it hangs for Grid 06.

I have used a bunch of exploratory cout statements, but I have not been able to figure out what is causing the program to hang.

Please let me know if you need more information to be able to understand my code.

#include <algorithm>
#include <array>
#include <fstream>
#include <iostream>
#include <map>
#include <set>
#include <vector>
using namespace std;

bool solve(map<int, set<int>>* X, map<int, array<int, 4>>* Y,
           vector<int>* solution);
void select_row(map<int, set<int>>* X, map<int, array<int, 4>>* Y,
                vector<set<int>>* cs, int r);
void deselect_row(map<int, set<int>>* X, map<int, array<int, 4>>* Y,
                  vector<set<int>>* cs, int r);

// square boxes only
static const int BOX_SIZE = 3; // standard sudoku grid
static const int SIZE = BOX_SIZE*BOX_SIZE;

int main() {
  int top_left_sum = 0;

  // initialize the sparse matrix representation of the sudoku problem
  map<int, set<int>> C; // constraints/columns
  for (int i = 0; i < 4*SIZE*SIZE; i++) {
    C[i] = set<int>();
  }
  map<int, array<int, 4>> R; // subsets/rows
  for (int i = 0; i < SIZE*SIZE*SIZE; i++) {
    // i is the subset index and encodes location and number on grid
    int index = i/SIZE;
    int row = i/(SIZE*SIZE);
    int column = (i/SIZE) % SIZE;
    int box = BOX_SIZE*(row/BOX_SIZE) + column/BOX_SIZE;
    int value = i % SIZE;

    // there are 4 constaints satisfied by each number placement
    array<int, 4> subset;  
    // insert the keys of constraints that subset satisfies
    subset[0] = (index); // row-column
    subset[1] = (SIZE*SIZE + SIZE*row + value); // row-number
    subset[2] = (2*SIZE*SIZE + SIZE*column + value); // column-number
    subset[3] = (3*SIZE*SIZE + SIZE*box + value); // box-number

    R.insert(pair<int, array<int, 4>>(i, subset));

    for (auto c : subset) {
      C[c].insert(i);
    }
  }

  ifstream ifs("../sudoku.txt");

  string line;
  while (getline(ifs, line)) {
    if (line[0] == 'G') {
      map<int, set<int>> X = C;
      map<int, array<int, 4>> Y = R;
      vector<int> solution;
      for (int i = 0; i < SIZE; i++) {
        getline(ifs, line);
        for (int j = 0; j < SIZE; j++) {
          if (line[j] != '0') {
            int r = SIZE*SIZE*i + SIZE*j + line[j] - '1';
            solution.push_back(r);
            vector<set<int>> cs;
            select_row(&X, &Y, &cs, r);
          }
        }
      }

      solve(&X, &Y, &solution);
      sort(solution.begin(), solution.end());

      top_left_sum += 100*(solution[0] % SIZE + 1) 
                      + 10*(solution[1] % SIZE + 1) 
                      + solution[2] % SIZE + 1;

      // display solution
      for (size_t i = 0; i < solution.size(); i++) {
        if (i % 9 == 0) cout << endl;
        cout << solution[i] % 9 + 1 << ' ';
      } cout << endl << endl;
    }
  }
  ifs.close();
  cout << top_left_sum << endl;
  return 0;
}

bool solve(map<int, set<int>>* X, map<int, array<int, 4>>* Y,
           vector<int>* solution) {
  if ((*X).empty()) return true;

  // find the column with the minimum number of nonzero elements
  map<int, set<int>>::iterator c_min = (*X).begin();
  for (map<int, set<int>>::iterator c = ++(*X).begin();
       c != (*X).end(); ++c) {
    if ((*c).second.size() < (*c_min).second.size()) {
      c_min = c;
    }
  }

  // for each row pointed to by c_min, call solve recursively
  for (auto r : (*c_min).second) {
    (*solution).push_back(r);
    vector<set<int>> cs;
    select_row(X, Y, &cs, r);
    if (solve(X, Y, solution)) return true;
    deselect_row(X, Y, &cs, r);
    (*solution).pop_back();
  }
  return false;
}

void select_row(map<int, set<int>>* X, map<int, array<int, 4>>* Y,
                vector<set<int>>* cs, int r) {
  for (auto j : (*Y)[r]) {
    for (auto i : (*X)[j]) {
      for (auto k : (*Y)[i]) {
        if (k != j) (*X)[k].erase(i);
      }
    }
    (*cs).push_back((*X)[j]);
    (*X).erase(j);
  }
  return;
}

void deselect_row(map<int, set<int>>* X, map<int, array<int, 4>>* Y,
                  vector<set<int>>* cs, int r) {
  for (array<int, 4>::reverse_iterator j = (*Y)[r].rbegin();
       j != (*Y)[r].rend(); ++j) {
    (*X)[*j] = (*cs).back();
    (*cs).pop_back();
    for (auto i : (*X)[*j]) {
      for (auto k : (*Y)[i]) {
        if (k != *j) (*X)[k].insert(i);
      }
    }
  }
  return;
}

Solution

As PaulMackenzie pointed out, the problem with my code was that I was erasing objects which still had pointers to them initialized, specifically the integers inside the set (*c_min).second that I iterate over in my solve function.

I fixed this by making a copy of (*c_min).second and iterating over the copy.

The fixed version of the solve function looks like this:

bool solve(map<int, set<int>>* X, map<int, array<int, 4>>* Y,
           vector<int>* solution) {
  if ((*X).empty()) return true;

  // find the column with the minimum number of nonzero elements
  map<int, set<int>>::iterator c_min = (*X).begin();
  for (map<int, set<int>>::iterator c = ++(*X).begin();
       c != (*X).end(); ++c) {
    if ((*c).second.size() < (*c_min).second.size()) {
      c_min = c;
    }
  }

  set<int> c = (*c_min).second; // ADDED LINE

  // for each row pointed to by c_min, call solve recursively
  for (auto r : c) { // CHANGED LINE
    (*solution).push_back(r);
    vector<set<int>> cs;
    select_row(X, Y, &cs, r);
    if (solve(X, Y, solution)) return true;
    deselect_row(X, Y, &cs, r);
    (*solution).pop_back();
  }
  return false;
}