Fast sudoku solver

solver.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <set>
#include <algorithm>

// first index is which row, second is which col;
#define GRID(I,J) grid[9*(I)+(J)]
#define GET(grid,I,J) (grid)[9*(I)+(J)]

std::set<int> all_numbers;

void zero_grid(int *grid){
	for(int i = 0; i < 81; ++i){
		grid[i] = 0;
	}
}

// returns how many already filled in
int read_initial_grid(int *grid, const char *filename){
	std::ifstream fin(filename);
	int count = 0;
	for(int i = 0; i < 9; ++i){
		for(int j = 0; j < 9; ++j){
			fin >> GRID(i,j);
			if(GRID(i,j) < 1 || GRID(i,j) > 9){
				GRID(i,j) = 0;
			}else{
				count++;
			}
		}
	}
	fin.close();
	return count;
}

void print_grid(int *grid){
	for(int i = 0; i < 9; ++i){
		for(int j = 0; j < 9; ++j){
			std::cout << " " << GRID(i,j);
		}
		std::cout << std::endl;
	}
}

void fill_grid(const char *reason, int *grid, int i, int j, int val){
	std::cout << "Filled in row " << i+1 << ", col " << j+1 << " = " << val << " by " << reason << std::endl;
	GRID(i,j) = val;
}

void generate_occlusion_map(int *grid, int *blocking_map, int n){
	// generated the grid of blocked cells (cells which cannot possibly contain value n
	for(int i = 0; i < 9; ++i){
		for(int j = 0; j < 9; ++j){
			// if the cell is already filled, then it is blocked
			if(GRID(i,j) != 0){
				GET(blocking_map,i,j) = 1;
			}
			
			if(GRID(i,j) == n){
				// if we found the number, then block out all rows/cols that contain it
				for(int k = 0; k < 9; ++k){
					GET(blocking_map,i,k) = 1;
					GET(blocking_map,k,j) = 1;
				}
				// also block out entire 3x3 block that contains it
				for(int k = 3*(i/3); k < 3*(i/3+1); ++k){
					for(int l = 3*(j/3); l < 3*(j/3+1); ++l){
						GET(blocking_map,k,l) = 1;
					}
				}
			}
		}
	}
}

int fill_by_block_using_occlusion_map(int *grid, int *blocking_map, int n, const char *reason){
	// Now see if any 3x3 blocks have only 1 zero
	for(int i = 0; i < 3; ++i){
		for(int j = 0; j < 3; ++j){ // these outer 2 loops iterate over each 3x3 block
		
			std::stringstream reason_str;
			reason_str << reason;
			
			int blocked_count = 0;
			int last_unblocked_k = -1, last_unblocked_l = -1;
			for(int k = 3*i+0; k < 3*i+3; ++k){
				for(int l = 3*j+0; l < 3*j+3; ++l){
					if(GET(blocking_map,k,l) == 1){
						blocked_count++;
					}else{
						last_unblocked_k = k;
						last_unblocked_l = l;
					}
				}
			}
			if(blocked_count == 8){
				// It's unique!
				reason_str << " in block " << i+1 << "," << j+1;
				fill_grid(reason_str.str().c_str(), grid, last_unblocked_k, last_unblocked_l, n);
				return 1;
			}
		}
	}
	return 0;
}

int fill_by_line_using_occlusion_map(int *grid, int *blocking_map, int n, const char *reason){
	for(int i = 0; i < 9; ++i){
		int num_holes_in_col_i = 0; // count number of unblocked cells (holes) in each col
		int num_holes_in_row_i = 0;
		int last_hole_in_col_i = -1;
		int last_hole_in_row_i = -1;
		for(int j = 0; j < 9; ++j){
			if(GET(blocking_map,i,j) == 0){
				num_holes_in_row_i++;
				last_hole_in_row_i = j;
			}
			if(GET(blocking_map,j,i) == 0){
				num_holes_in_col_i++;
				last_hole_in_col_i = j;
			}
		}
		std::stringstream reason_str;
		reason_str << reason;
		if(num_holes_in_row_i == 1){
			reason_str << " on row " << i+1;
			fill_grid(reason_str.str().c_str(), grid, i, last_hole_in_row_i, n);
			return 1;
		}
		if(num_holes_in_col_i == 1){
			reason_str << " on col " << i+1;
			fill_grid(reason_str.str().c_str(), grid, last_hole_in_col_i, i, n);
			return 1;
		}
	}
	return 0;
}

// First order occlusion test
// returns how many could be filled in
int occlusion_test(int *grid, int n){
	int b[81]; zero_grid(b);
	int ret = 0;
	
	generate_occlusion_map(grid, b, n);
	
	ret += fill_by_block_using_occlusion_map(grid, b, n, "occlusion test");
	if(ret == 0){
		ret += fill_by_line_using_occlusion_map(grid, b, n, "occlusion test");
	}
	return ret;
}

int block_refine_occlusion_map(int *b){
	int refinements = 0;
	// Now see if any 3x3 blocks have unblocked cells in only a single row/col
	for(int i = 0; i < 3; ++i){
		for(int j = 0; j < 3; ++j){ // these outer 2 loops iterate over each 3x3 block
			int min_unblocked_k = 10;
			int max_unblocked_k = -1;
			int min_unblocked_l = 10;
			int max_unblocked_l = -1;
			for(int k = 3*i+0; k < 3*i+3; ++k){
				for(int l = 3*j+0; l < 3*j+3; ++l){
					if(GET(b,k,l) == 0){
						if(k < min_unblocked_k){ min_unblocked_k = k; }
						if(k > max_unblocked_k){ max_unblocked_k = k; }
						if(l < min_unblocked_l){ min_unblocked_l = l; }
						if(l > max_unblocked_l){ max_unblocked_l = l; }
					}
				}
			}
			
			if(min_unblocked_k == max_unblocked_k){ // all in single row
				for(int k = 0; k < 9; ++k){
					if(3*j+0 <= k && k < 3*j+3){ continue; } // remember to only add occlusion in OTHER blocks
					if(GET(b,min_unblocked_k,k) == 0){
						++refinements;
					}
					GET(b,min_unblocked_k,k) = 1;
				}
			}else if(min_unblocked_l == max_unblocked_l){
				for(int k = 0; k < 9; ++k){
					if(3*i+0 <= k && k < 3*i+3){ continue; } // remember to only add occlusion in OTHER blocks
					if(GET(b,k,min_unblocked_l) == 0){
						++refinements;
					}
					GET(b,k,min_unblocked_l) = 1;
				}
			}
		}
	}
	return refinements;
}

// Second order occlusion test
// Perform first order occlusion, then if any 3x3 blocks contain unblocked cells in a single row/col, use that to perform another level of occlusion
int occlusion_test2(int *grid, int n){
	int b[81]; zero_grid(b);
	int ret = 0;
	
	generate_occlusion_map(grid, b, n);

	// It is assumed that a first order occlusion test was performed already
	
	// Now see if any 3x3 blocks have unblocked cells in only a single row/col
	int p = 1;
	while(block_refine_occlusion_map(b)){
		++p;
	}
	
	std::stringstream reason;
	reason << "order " << p << " occlusion test";
	
	ret += fill_by_block_using_occlusion_map(grid, b, n, reason.str().c_str());
	if(ret == 0){
		ret += fill_by_line_using_occlusion_map(grid, b, n, reason.str().c_str());
	}
	return ret;
}

// Use standard occlusion test, but also go through all numbers != n to see if their feasible locations match up exactly, then those numbers must be confined to the feasible locations, which act as a group to block n from being in those locations
int occlusion_test_with_feasible_set_reduction(int *grid, int n){
	int b[81]; zero_grid(b);
	int ret = 0;
	
	generate_occlusion_map(grid, b, n);

	// generate occlusion maps for all other numbers
	int b2[9][81];
	for(int n2 = 1; n2 <= 9; ++n2){ if(n != n2){
		zero_grid(b2[n2-1]);
		generate_occlusion_map(grid, b2[n2-1], n2);
		
	}}
	
	std::stringstream sout;
	sout << "occlusion test by feasible set reduction, sets: ";
	
	for(int i = 0; i < 3; ++i){
		for(int j = 0; j < 3; ++j){ // these outer 2 loops iterate over each 3x3 block
			// find which numbers are already in this block, so we can skip it for all n2 already in this block
			int numbers_in_block[9];
			int num_numbers_in_block = 0;
			for(int k = 3*i+0; k < 3*i+3; ++k){
				for(int l = 3*j+0; l < 3*j+3; ++l){
					if(GRID(k,l) != 0){
						numbers_in_block[num_numbers_in_block] = GRID(k,l);
						num_numbers_in_block++;
					}
				}
			}
		
			// See if occlusion maps of other numbers line up exactly within any 3x3 block
			// easiest way is to xor the maps; if it comes up all zeros within a block then they are identical
			// these xors must be performed pairwise, so we get a (strictly upper triangular) matrix of binary matched or not values
			int matches[81]; zero_grid(matches);
			int feasible_set_size[9]; // number of 0's for each n2
			// We have a feasible set overlap of size 2 if matches(i,j) = 1, feasible_set_size(i)=feasible_set_size(j) = 2
			// We have a feasible set overlap of size 3 if matches(i,j) = matches(j,k) = 1, feasible_set_size(i) = feasible_set_size(j) = feasible_set_size(k) = 3
			// etc.
			
			// fill in feasible_set_size
			for(int n2 = 1; n2 <= 9; ++n2){
				bool already_in = false;
				for(int nn = 0; nn < num_numbers_in_block; ++nn){ if(n2 == numbers_in_block[nn]){ already_in = true; }}
				if(n != n2 && !already_in){
					feasible_set_size[n2-1] = 0;
					for(int k = 3*i+0; k < 3*i+3; ++k){
						for(int l = 3*j+0; l < 3*j+3; ++l){
							if(GET(b2[n2-1],k,l) == 0){
								feasible_set_size[n2-1]++;
							}
						}
					}
				}
			}
			
			// fill in matches
			for(int n2 = 1; n2 <= 9; ++n2){
				bool already_in = false;
				for(int nn = 0; nn < num_numbers_in_block; ++nn){ if(n2 == numbers_in_block[nn]){ already_in = true; }}
				if(n != n2 && !already_in){
				
					for(int n3 = n2+1; n3 <= 9; ++n3){
						bool already_in2 = false;
						for(int nn2 = 0; nn2 < num_numbers_in_block; ++nn2){ if(n3 == numbers_in_block[nn2]){ already_in2 = true; }}
						if(n != n3 && !already_in2){
							// perform the xor; count number of differences instead
							int num_differences = 0;
							for(int k = 3*i+0; k < 3*i+3; ++k){
								for(int l = 3*j+0; l < 3*j+3; ++l){
									if(GET(b2[n2-1],k,l) != GET(b2[n3-1],k,l)){
										num_differences++;
										k = 9; break; // exit loop with variable k
									}
								}
							}

							// At this point, if num_differences > 0, the XOR was not all zero
							if(num_differences > 0){
								GET(matches,n2-1,n3-1) = 0;
							}else{
								GET(matches,n2-1,n3-1) = 1;
							}
						}
					}
				}
			}
			
			// Now test to see if there are collections of feasible sets that did overlap
			// count how many on each row of matches contains a 1
			for(int row = 0; row < 9; ++row){
				bool already_in = false;
				for(int nn = 0; nn < num_numbers_in_block; ++nn){ if(row+1 == numbers_in_block[nn]){ already_in = true; }}
				if(n != row+1 && !already_in){
					
					int num_overlaps = 1; // at least occlusion map of row+1 overlaps with itself
					int overlapping_numbers[9]; // keep track of which numbers overlapped
					overlapping_numbers[0] = row+1;
					for(int col = row+1; col < 9; ++col){
						if(GET(matches,row,col) == 1){
							overlapping_numbers[num_overlaps] = col+1;
							num_overlaps++;
						}
					}
					// we limit ourselves to sizes of 2 through 4 to be reasonable
					if(2 <= num_overlaps && num_overlaps <= 4){ // subtract 1 since the 1's on the diagonals are never counted
						bool can_reduce = true;
						for(int n_overlaps = 0; n_overlaps < num_overlaps; ++n_overlaps){
							if(feasible_set_size[overlapping_numbers[n_overlaps]-1] != num_overlaps){
								can_reduce = false;
								break;
							}
						}
						if(can_reduce){
							// now take the feasible set (1-b2[row][]) and count it as blocking
							for(int k = 3*i+0; k < 3*i+3; ++k){
								for(int l = 3*j+0; l < 3*j+3; ++l){
									if(GET(b2[row], k, l) == 0){
										GET(b, k, l) = 1;
									}
								}
							}
							sout << "(" << overlapping_numbers[0];
							for(int n_overlaps = 1; n_overlaps < num_overlaps; ++n_overlaps){
								sout << "," << overlapping_numbers[n_overlaps];
							}
							sout << ") in block " << i+1 << "," << j+1 << " ";
						}
					}
				}
			}
		}
	}
	
	ret += fill_by_block_using_occlusion_map(grid, b, n, sout.str().c_str());
	if(ret == 0){
		ret += fill_by_line_using_occlusion_map(grid, b, n, sout.str().c_str());
	}
	return ret;
}

int trivial_completion(int *grid){
	// go through each block
	for(int i = 0; i < 3; ++i){
		for(int j = 0; j < 3; ++j){ // these outer 2 loops iterate over each 3x3 block
			int block_numbers[9];
			int num_block_numbers = 0;
			int last_hole_k = -1;
			int last_hole_l = -1;
			for(int k = 3*i+0; k < 3*i+3; ++k){
				for(int l = 3*j+0; l < 3*j+3; ++l){
					if(GRID(k,l) != 0){
						block_numbers[num_block_numbers] = GRID(k,l);
						num_block_numbers++;
					}else{
						last_hole_k = k;
						last_hole_l = l;
					}
				}
			}
			if(num_block_numbers == 8){
				for(int n = 1; n <= 9; ++n){
					bool in_block = false;
					for(int k = 0; k < num_block_numbers; ++k){
						if(block_numbers[k] == n){
							in_block = true;
							break;
						}
					}
					if(!in_block){
						std::stringstream reason_str;
						reason_str << "trivial completion in block " << i+1 << "," << j+1;
						fill_grid(reason_str.str().c_str(), grid, last_hole_k, last_hole_l, n);
						return 1;
					}
				}
			}
		}
	}
	
	// go through each line
	for(int i = 0; i < 9; ++i){
		int num_holes_in_row_i = 0;
		int num_holes_in_col_i = 0;
		int last_hole_j_in_row_i = -1;
		int last_hole_j_in_col_i = -1;
		for(int j = 0; j < 9; ++j){
			if(GRID(i,j) == 0){
				num_holes_in_row_i++;
				last_hole_j_in_row_i = j;
			}
			if(GRID(j,i) == 0){
				num_holes_in_col_i++;
				last_hole_j_in_col_i = j;
			}
		}
		if(num_holes_in_row_i == 1){
			// find which number it was
			for(int n = 1; n <= 9; ++n){
				bool in_line = false;
				for(int j = 0; j < 9; ++j){
					if(GRID(i,j) == n){
						in_line = true;
						break;
					}
				}
				if(!in_line){
					std::stringstream reason_str;
					reason_str << "trivial completion in row " << i+1;
					fill_grid(reason_str.str().c_str(), grid, i, last_hole_j_in_row_i, n);
					return 1;
				}
			}
		}else if(num_holes_in_col_i == 1){
			// find which number it was
			for(int n = 1; n <= 9; ++n){
				bool in_line = false;
				for(int j = 0; j < 9; ++j){
					if(GRID(j,i) == n){
						in_line = true;
						break;
					}
				}
				if(!in_line){
					std::stringstream reason_str;
					reason_str << "trivial completion in col " << i+1;
					fill_grid(reason_str.str().c_str(), grid, last_hole_j_in_col_i, i, n);
					return 1;
				}
			}
		}
	}
	return 0;
}

// go through each empty cell and find which numbers cannot go in that cell
// any numbers in the row, column, or block are eliminated
int cell_wise_elimination(int *grid){
	for(int i = 0; i < 9; ++i){
		for(int j = 0; j < 9; ++j){
			if(GRID(i,j) == 0){
				// NOw, we have found an empty cell
				std::set<int> eliminated_numbers;
				// check block
				for(int k = 3*(i/3); k < 3*(i/3+1); ++k){
					for(int l = 3*(j/3); l < 3*(j/3+1); ++l){
						if(GRID(k,l) != 0){
							eliminated_numbers.insert(GRID(k,l));
						}
					}
				}
				// check row/col
				for(int k = 0; k < 9; ++k){
					if(GRID(i,k) != 0){
						eliminated_numbers.insert(GRID(i,k));
					}
					if(GRID(k,j) != 0){
						eliminated_numbers.insert(GRID(k,j));
					}
				}
				if(eliminated_numbers.size() == 8){
					std::set<int> C;
					std::set_difference(all_numbers.begin(), all_numbers.end(), eliminated_numbers.begin(), eliminated_numbers.end(), std::inserter(C, C.begin()));
					int n = *(C.begin());
					fill_grid("cell-wise elimination", grid, i, j, n);
					return 1;
				}
			}
		}
	}
	return 0;
}

// look at row,col,block and find which numbers can go in which cells
int feasible_set_analysis(int *grid){
	int bn[9][81];
	for(int n = 1; n <= 9; ++n){
		zero_grid(bn[n-1]);
		
		// Generate order-p occlusion map
		generate_occlusion_map(grid, bn[n-1], n);
		int p = 1;
		while(block_refine_occlusion_map(bn[n-1])){
		}
	}
	
	std::set<int> feas[81];
	int max_set_size = 3;
	for(int i = 0; i < 9; ++i){
		for(int j = 0; j < 9; ++j){
			if(GRID(i,j) == 0){
				for(int n = 1; n <= 9; ++n){
					if(GET(bn[n-1],i,j) == 0){
						GET(feas,i,j).insert(n);
					}
				}
				if(max_set_size < GET(feas,i,j).size()){ max_set_size = GET(feas,i,j).size(); }
			}
		}
	}
	
	// print out feasible sets
	for(int i = 0; i < 9; ++i){
		for(int j = 0; j < 9; ++j){
			if(GRID(i,j) == 0){
				int num_printed = 0;
				for(std::set<int>::const_iterator it = GET(feas,i,j).begin(); it != GET(feas,i,j).end(); ++it){
					std::cout << *it;
					++num_printed;
				}
				for(int k = num_printed; k <= max_set_size; ++k){
					std::cout << " ";
				}
			}else{
				std::cout << "[" << GRID(i,j) << "]";
				for(int k = 3; k <= max_set_size; ++k){
					std::cout << " ";
				}
			}
		}
		std::cout << std::endl;
	}
}


#include "brute.cpp"

int main(int argc, char *argv[]){
	int grid[81];
	
	if(argc < 2){
		std::cout << "Usage: " << argv[0] << " input.txt" << std::endl;
		return 0;
	}
	
	int initial_count = read_initial_grid(grid, argv[1]);
	int total_count = initial_count;

	for(int n = 1; n <= 9; ++n){
		all_numbers.insert(n);
	}
	
	int num_filled_per_pass;
	do{
		num_filled_per_pass = 0;
		if(num_filled_per_pass == 0){
			num_filled_per_pass += trivial_completion(grid);
		}
		for(int n = 1; n <= 9; ++n){ if(num_filled_per_pass == 0){
			num_filled_per_pass += occlusion_test(grid, n);
		}}
		for(int n = 1; n <= 9; ++n){ if(num_filled_per_pass == 0){
			num_filled_per_pass += occlusion_test2(grid, n);
		}}
		for(int n = 1; n <= 9; ++n){ if(num_filled_per_pass == 0){
			num_filled_per_pass += occlusion_test_with_feasible_set_reduction(grid, n);
		}}
		if(num_filled_per_pass == 0){
			num_filled_per_pass += cell_wise_elimination(grid);
		}
		total_count += num_filled_per_pass;
	}while(num_filled_per_pass > 0 && total_count < 81);
	
	std::cout << "Total filled in: " << total_count << std::endl;
	
	print_grid(grid);
	
	if(total_count < 81){
		std::cout << "Built-in rules could not complete the solve. Remaining feasible sets:" << std::endl;
		feasible_set_analysis(grid);
		std::cout << "Transfering to brute-force solver" << std::endl;
		
		Puzzle p;
		for(int i = 0; i < 9; i++){
			for(int j = 0; j < 9; j++){
				if(GRID(i,j) != 0){
					p.set(i,j,GRID(i,j));
				}
			}
		}
		if(solve(p, 0, gridSize, 0, gridSize, true, true)){
			printPuzzle(p);
		}else{
			std::cout << "Failed to solve" << std::endl;
		}
	}
	
	return 0;
}