xsudoku.cpp

数独并行算法 intel 多线程优化大赛参赛作品

/*/////////////////////////////////////////
// 版权所有(C)  2000-2008 邓辉           //
// Email：      [e-mail]  //
// 说明：       Intel线程优化大赛参赛作品//
/////////////////////////////////////////*/

#include "XSudoku.h"

// 包含1的个数 
// 小于XFILLEDMASK为正常值
// 大于等于XFILLEDMASK值均为255  [XFindMinFeasibilityPos中减少判断]
__declspec(align(16)) uint8  g_xPopCnt[0xFFFF] = {0};
// 单元格的值到输出字符的映射
__declspec(align(16)) char   g_xOutTab[1 << XGRID_BSIZE] = {0};
// 输出缓冲
__declspec(align(16)) char   g_xOutCells[XGRID_BSIZE + XGRID_BW + 1][(XGRID_BSIZE + XGRID_BH) * 2 + 2] = {0};

// 保存与当前单元格相关的单元格
__declspec(align(16)) uint8  g_xMutuality[XGRID_AREA][XGRID_MUTUALITYSSE] = {0};
__declspec(align(16)) uint8  g_nMutuality[XGRID_AREA] = {0};

// 初始化静态加速Tab
void XInitTab()
{
	int i,j;

	// 初始化PopCnt数组
	g_xPopCnt[0] = 0;
	g_xPopCnt[1] = 1;
	g_xPopCnt[2] = 1;
	g_xPopCnt[3] = 2;

	for(i = 0x04;  i < 0x10;  ++i)   g_xPopCnt[i] = g_xPopCnt[i >> 2] + g_xPopCnt[i & 0x03];
	for(i = 0x10;  i < 0x100; ++i)   g_xPopCnt[i] = g_xPopCnt[i >> 4] + g_xPopCnt[i & 0x0F];
#pragma omp parallel for
	for(i = 0x100; i < XFILLEDMASK; ++i) g_xPopCnt[i] = g_xPopCnt[i >> 8] + g_xPopCnt[i & 0xFF];
#pragma omp parallel for
	for(i = XFILLEDMASK; i < 0x10000; ++i) g_xPopCnt[i] = 0xFF;

	// 计算每个单元格的关联单元格数组
	for(i = 0; i < XGRID_AREA; ++i)
	{
		g_nMutuality[i] = 0;

		for(j = 0; j < XGRID_AREA; ++j)
		{
			if(i == j) continue;

			// 如果两个单元格满足以下条件 则这两个单元格存在关联关系
			// 1. 单元格的横坐标相同 
			// 2. 单元格的纵坐标相同
			// 3. 单元格所在块的横坐标相同且单元格所在块的纵坐标相同[两个单元格在同一个块中]
			if( ( i % XGRID_BSIZE == j % XGRID_BSIZE ) || ( i / XGRID_BSIZE == j / XGRID_BSIZE ) || 
				( ( i / XGRID_BSIZE / XGRID_BH  == j / XGRID_BSIZE / XGRID_BH ) && 
				( ( i % XGRID_BSIZE ) / XGRID_BW ) == ( ( j % XGRID_BSIZE ) / XGRID_BW ) ) )
			{
				g_xMutuality[i][g_nMutuality[i]++] = j;
			}
		}

		for(j = XGRID_MUTUALITY; j < XGRID_MUTUALITYSSE; ++j)  g_xMutuality[i][j] = XGRID_AREA;
	}

	// 初始化值到字符的映射
	for(i = 0; i < XGRID_BSIZE; ++i) g_xOutTab[ 1 << i] = i + '1';

	// 初始化输出字符表格
	for(i = 0; i < XGRID_BSIZE + XGRID_BW + 1; ++i)
	{
		if(i % (XGRID_BH + 1) == 0)
		{
			for(j = 0; j < (XGRID_BSIZE + XGRID_BH) * 2 + 1; ++j)
			{
				g_xOutCells[i][j] =  (j % (XGRID_BW * 2 + 2) == 0) ? '+' : '-';
			}
		}
		else
		{
			for(j = 0; j < (XGRID_BSIZE + XGRID_BH) * 2 + 1; ++j)
			{
				g_xOutCells[i][j] =  (j % (XGRID_BW * 2 + 2) == 0) ? '|' : ' ';
			}
		}

		g_xOutCells[i][(XGRID_BSIZE + XGRID_BH) * 2 + 1] =  (i == XGRID_BSIZE + XGRID_BW) ? '\0' : '\n';
	}
}

// 输出网格
__inline void XAddSolution(int* pRunMode, int* pTotal, XCells pCells)
{
#ifndef XOUT_SOLUTION
	if(++(*pTotal) > 1)
	{
		*pRunMode &= ~XRUN_FIND_FLAG;
	}
#else
	static mutex xMutex;

	if(*pRunMode & XRUN_FIND_ONE)
	{
		if(++(*pTotal) > 1)
		{
			*pRunMode &= ~XRUN_FIND_FLAG;
		}
	}
	else
	{
		tbb::mutex::scoped_lock lock(xMutex);

		++(*pTotal);
	}

	if(*pRunMode & (XRUN_OUT_STD | XRUN_OUT_FILE))
	{
		for(int nPos = 0; nPos < XGRID_AREA; ++nPos)
		{
			int x = nPos % XGRID_BSIZE;
			int y = nPos / XGRID_BSIZE;
			g_xOutCells[y + y / XGRID_BH + 1][(x + x / XGRID_BW + 1) * 2] = g_xOutTab[pCells[nPos] - XFILLEDMASK];
		}

		if(*pRunMode & XRUN_OUT_STD) 
		{
			printf("第%d种\n%s\n\n", *pTotal, (const char*)g_xOutCells);
		}
		else
		{
			FILE* fp = fopen("xout.txt", "a");
			fprintf(fp, "第%d种\n%s\n\n", *pTotal, (const char*)g_xOutCells);
			fclose(fp);
		}
	}
#endif
}


// 为nPos位置的单元格置上nMask表示的值
// 当某个关联单元格清除后仅剩下一种可能性则继续递归调用XFillCell填充掉该单元格
// 将填充的单元格位置保存到pGrid->pFilledPos中
__inline int XFillCell(XGrid* pGrid, XCells pCells, int nPos, int nMask)
{
	int nFilled = pGrid->nFilled;

	for(;;)
	{
#ifdef XGRID_MUTUALITYSSE_CLEAR
		uint8* const pMutuality = pGrid->pMutuality[nPos];
		const uint8  nMutuality = pGrid->nMutuality[nPos];
#else
		uint8* const pMutuality = g_xMutuality[nPos];
		const uint8  nMutuality = g_nMutuality[nPos];
#endif

		pCells[nPos] = nMask + XFILLEDMASK;

		pGrid->pFilledPos[nFilled++] = nPos;

		nPos = -1;

#pragma unroll(16)
		for(int i = 0; i < nMutuality; ++i)
		{
			const uint8 nTmp = pMutuality[i];

			// 去掉该单元格填充nMask的可能性
			pCells[nTmp] &= ~nMask;

			// 如果单元格没有可选择的填值方案 返回失败
			if(pCells[nTmp] == 0) return XFILL_INVALID;

			// 判断该单元格是否只有一种填值可能,如果是将调用XFillCell给该单元格填值
			if(g_xPopCnt[pCells[nTmp]] == 1) nPos = nTmp;
		}

		if(nPos < 0) break;

		nMask = pCells[nPos];
	}

	pGrid->nFilled = nFilled;
	
	return pGrid->nFilled;
}

__inline int XFillCellLoad(XGrid* pGrid, XCells pCells, int nPos, int nMask)
{
	int nNewPos = -1;
#ifdef XGRID_MUTUALITYSSE_CLEAR
	uint8* const pMutuality = pGrid->pMutuality[nPos];
	const uint8  nMutuality = pGrid->nMutuality[nPos];
#else
	uint8* const pMutuality = g_xMutuality[nPos];
	const uint8  nMutuality = g_nMutuality[nPos];
#endif

	// 如果该单元格已经填充了nMask 直接返回
	if(pCells[nPos] == nMask + XFILLEDMASK) return pGrid->nFilled;

	// 判断当前单元格能填nMask表示的值
	if((pCells[nPos] & nMask) == 0) return XFILL_INVALID;

	// 为第nPos个单元格填上nMask表示的值
	pCells[nPos] = nMask + XFILLEDMASK;

	for(int i = 0; i < nMutuality; ++i)
	{
		const uint8 nTmp = pMutuality[i];

		// 去掉该单元格填充nMask的可能性
		pCells[nTmp] &= ~nMask;

		// 如果单元格没有可选择的填值方案 返回失败
		if(pCells[nTmp] == 0) return XFILL_INVALID;

		// 判断该单元格是否只有一种填值可能,如果是将调用XFillCell给该单元格填值
		if(g_xPopCnt[pCells[nTmp]] == 1) nNewPos = nTmp;
	}

	// 记录当前填充的单元格
	pGrid->pFilledPos[pGrid->nFilled++] = nPos;

	if(nNewPos >= 0) 
	{
		if(XFillCell(pGrid, pCells, nNewPos, pCells[nNewPos]) == XFILL_INVALID)
		{
			// 恢复已经填值的单元格
			--pGrid->nFilled;

			return XFILL_INVALID;
		}
	}

	return pGrid->nFilled;
}


// 查找可填的数字最少的单元格
__inline int XFindMinFeasibilityPos(XCells pCells)
{
	int nMinCnt = 0xFF;
	int nMinPos = -1;

	// 在所有未填值的单元格中查找可填值方案最少的单元格
	for(int i = 0; i < XGRID_AREA && nMinCnt > 1; ++i)
	{
		if(g_xPopCnt[pCells[i]] < nMinCnt)
		{
			nMinPos = i;
			nMinCnt = g_xPopCnt[pCells[i]];
		}
	}

	return nMinPos;
}

// 删除与单元格nPos相关联的单元格 [XGRID_BSIZE>=9且需要查找所有情况时清除关联可加速]
#ifdef XGRID_MUTUALITYSSE_CLEAR
	__inline void XRemoveMutuality(XGrid* pGrid, uint8 nPos)
	{
		uint8* const pCur = pGrid->pMutuality[nPos];
		for(int i = pGrid->nMutuality[nPos] - 1; i >= 0; --i)
		{
			uint8* const pTmp = pGrid->pMutuality[pCur[i]];
			const int  nTmp = --pGrid->nMutuality[pCur[i]];
			for(int j = 0; j < nTmp; ++j)
			{
				if(pTmp[j] == nPos)
				{
					for(int k = j; k < nTmp; ++k) pTmp[k] = pTmp[k+1];
					break;
				}
			}
		}
	}
#endif

// 加载测试数据
__inline int XLoadGrid(XGrid* pGrid, const char* pInput)
{
	int nPos = 0;

	int nMask = (1 << XGRID_BSIZE) - 1;

#ifdef XGRID_MUTUALITYSSE_CLEAR
	XMEMCOPY(pGrid->pMutuality, g_xMutuality, sizeof(g_xMutuality));
	XMEMCOPY(pGrid->nMutuality, g_nMutuality, sizeof(g_nMutuality));
#endif

	for(nPos = 0; nPos < XGRID_AREA; ++nPos) pGrid->xCells[nPos] = nMask;

	for(nPos = 0; nPos < XGRID_AREA; ++nPos)
	{
		if(pInput[nPos] >= '1' && pInput[nPos] <= '0' + XGRID_BSIZE)
		{	// 填充nPos的值为 1 << ( pInput[nPos] - '1' )
			if(XFillCellLoad(pGrid, pGrid->xCells, nPos, 1 << ( pInput[nPos] - '1' )) == XFILL_INVALID)
			{
				return XFILL_INVALID;
			} 
		}
	}

#ifdef XGRID_MUTUALITYSSE_CLEAR
	// [XGRID_BSIZE>=9且需要查找所有情况时清除关联可加速]
	for(int i = 0; i < pGrid->nFilled; ++i)
	{
		XRemoveMutuality(pGrid, pGrid->pFilledPos[i]);
	}
#endif

	return XGRID_AREA - pGrid->nFilled;
}
 
__inline void XFillCellStack(XGrid* pGrid)
{
	__declspec(align(16)) XCells    xCellsStack[XGRID_AREA  + 1];
	__declspec(align(16)) int		nMinPosStack[XGRID_AREA];
	__declspec(align(16)) int		nCurValStack[XGRID_AREA];
	__declspec(align(16)) int		nFilledStack[XGRID_AREA];
	int nIndex = 0;

	XMEMCOPY(xCellsStack[0], pGrid->xCells, sizeof(xCellsStack[0]));

	// 查找当前填值可能性最少的单元格
	nMinPosStack[nIndex] = XFindMinFeasibilityPos(xCellsStack[nIndex]);
	nCurValStack[nIndex] = xCellsStack[nIndex][nMinPosStack[nIndex]];
	nFilledStack[nIndex] = pGrid->nFilled;

	while( (*pGrid->pRunMode & XRUN_FIND_FLAG) && nIndex >= 0 )
	{
		nFilledStack[nIndex + 1] = XFILL_INVALID;

		// 尝试填充当前单元格
		while(nFilledStack[nIndex + 1] == XFILL_INVALID && nCurValStack[nIndex])
		{
			XMEMCOPY(xCellsStack[nIndex + 1], xCellsStack[nIndex], sizeof(xCellsStack[0]));
			unsigned int nMask = nCurValStack[nIndex] & -nCurValStack[nIndex];
			nCurValStack[nIndex] -= nMask;
			nFilledStack[nIndex + 1] = XFillCell(pGrid, xCellsStack[nIndex + 1], nMinPosStack[nIndex], nMask);
		}

		if(nFilledStack[nIndex + 1] > XGRID_AREA)
		{	// 没有合法的填充值
			pGrid->nFilled = nFilledStack[--nIndex];
		}
		else if(nFilledStack[nIndex + 1] != XGRID_AREA)
		{	// 是合法的填充值, 还存在未填值的单元格
			++nIndex;
			// 查找当前填值可能性最少的单元格
			nMinPosStack[nIndex] = XFindMinFeasibilityPos(xCellsStack[nIndex]);
			nCurValStack[nIndex] = xCellsStack[nIndex][nMinPosStack[nIndex]];
		}
		else
		{	// 所有单元格均已填值
			XAddSolution(pGrid->pRunMode, pGrid->pTotal, xCellsStack[nIndex + 1]);
			pGrid->nFilled = nFilledStack[nIndex];
		}
	}
}

class CXSudokuTask: public tbb::task 
{
	XGrid*	m_pGrid;
public:
	CXSudokuTask( XGrid* pGrid) : m_pGrid(pGrid){}

	tbb::task* execute() 
	{
		if( ( *m_pGrid->pRunMode & XRUN_FIND_FLAG ) == 0) return NULL;

		if( m_pGrid->nFilled < XGRID_CHANGE_MUTUALITY_SIZE) 
		{
			tbb::task_list list;

			__declspec(align(16)) XGrid  pChildGrid[XGRID_BSIZE];

			int nChildGrid = 0, nNewFilled = 0, nOldFilled = m_pGrid->nFilled;

			// 查找当前填值可能性最少的单元格
			int nMinPos = XFindMinFeasibilityPos(m_pGrid->xCells);

			// 如果当前单元格只有一种填值情况,则直接给该单元格填值
			if(g_xPopCnt[m_pGrid->xCells[nMinPos]] == 1 && nOldFilled < XGRID_AREA)
			{
				nNewFilled = XFillCell(m_pGrid, m_pGrid->xCells, nMinPos, m_pGrid->xCells[nMinPos]);
				nMinPos = XFindMinFeasibilityPos(m_pGrid->xCells);
			}

			if(nNewFilled < XGRID_AREA)
			{	// 存在多种可能的填值情况
				int  nCurVal = m_pGrid->xCells[nMinPos];

				while(nCurVal)
				{
					XGrid* pCurGrid = pChildGrid + nChildGrid;

					XMEMCOPY(pCurGrid, m_pGrid, sizeof(pChildGrid[0]));

					unsigned int nMask = nCurVal & -nCurVal;

					nCurVal -= nMask;

					// 给该单元格填可能的值
					nNewFilled = XFillCell(pCurGrid, pCurGrid->xCells, nMinPos, nMask);

					if(nNewFilled < XGRID_AREA) 
					{	// 是合法的填充值, 还存在未填值的单元格
						#ifdef XGRID_MUTUALITYSSE_CLEAR
							// [XGRID_BSIZE>=9且需要查找所有情况时清除关联可加速]
							for(int i = nOldFilled; i < nNewFilled; ++i)
							{
								XRemoveMutuality(pCurGrid, pCurGrid->pFilledPos[i]);
							}
						#endif

						// 增加子Task
						list.push_back( *new( allocate_child() ) CXSudokuTask(pCurGrid));

						++nChildGrid;
					}
					else if(nNewFilled == XGRID_AREA)
					{	// 所有单元格均已填值
						XAddSolution(pCurGrid->pRunMode, pCurGrid->pTotal, pCurGrid->xCells);
					}
				}

				set_ref_count( nChildGrid + 1 );

				spawn_and_wait_for_all(list);
			}
			else if(nNewFilled == XGRID_AREA)
			{	// 所有单元格均已填值
				XAddSolution(m_pGrid->pRunMode, m_pGrid->pTotal, m_pGrid->xCells);
			}
		} 
		else
		{	
			XFillCellStack(m_pGrid); 
		}

		return NULL;
	}
};

// 大于 9 * 9 的网格 查找所有的填值可能方案时 效率更高 
void XParallelSudoku( XGrid* pGrid ) 
{
	CXSudokuTask& xTask = *new(tbb::task::allocate_root()) CXSudokuTask(pGrid);

	tbb::task::spawn_root_and_wait(xTask);
}

int	XFindGrid(int nRunMode, const char* pInput)
{
	int nTotal = 0;

	XGrid xGrid = {0};

	xGrid.pRunMode = &nRunMode;
	xGrid.pTotal   = &nTotal;

	// 加载网格数据
	switch(XLoadGrid(&xGrid, pInput))
	{
	case 0:
		nTotal = 1;
		break;
	case XFILL_INVALID:
		break;
	default:
		// 填充单元格
		XFillCellStack(&xGrid);
	}

	return nTotal;
}