对于数组使用外存的探索

最近遇到一个题，题目是这样描述的：对于一个数n，如果n是偶数，则导出n / 2，否则导出3 * n + 1，依次类推，直到最后结果是1，由此得到一个以n起始的Collatz数列。归纳为数学公式如下：

n → n / 2, n is even
n → 3 * n + 1, n is odd

举个例子，5→16→8→4→2→1，则以5为起始的Collatz数列长度为6。
当然题目难在数很大：要找出以1000000以内的数为开头的最长的Collatz数列的长度。

算法与数据结构

很明显，100万不小，如果单穷举的话时间肯定不能接受。还因为数链有交叉，所以我想到了动态规划，还因为中间算数过程会出现大大超出100万的现象，我决定数组使用外存。为了方便演示，暂时只考虑1000以内的数。我是具体如下设计的：

const int LEN = 1000;
const int SIZE = 4 * LEN; 
vector < int >cnt(SIZE);

这样，cnt就存储了SIZE个数的Collatz 数组长度。SIZE是LEN的4倍也仅仅是因为它能保证cnt[LEN-1]位导出cnt[3(LEN-1) +1] 不越界，这里仅作参考。
接下来就是外存的设计，我把数组以二进制形式存储在文件中，每次需要特定范围的数时就读取出相应的文件。文件命名也有规律，比如文件都存在文件夹./Euler14-data/中，文件名前缀统一为Euler14-，后缀统一为.data，文件名格式就如Euler14-层号.data，层号从0开始，第0层就包括0～SIZE-1，第1层就包括SIZE～2SIZE-1，依次类推。换句话说是cnt[i]表示的其实是SIZE*层号+i的Collatz数列长度。另外我定义了一个全局变量totalSize用来记录这些文件一共的存储空间大小。
然后就是对文件的读取和写入了，我将它们写成void download(const string & name)与void storage(const string & name)函数了。只需要把需要的对应层号的文件名传给&name即可。
紧接着就是核心算法的设计：

int collatz(const int layer, const int i, const int preview)	/*layer表示当前层，preview表示上一层，i表示当前层中位置*/
{
…
        /*判断是否是同层，如果不是，就要存储上一层，下载当前层了*/
	if (layer != preview)
	{
		if (preview >= 0)
		    storage(preName);
		download(curName);
	}
        /*边界条件*/
	if (layer == 0 && (i == 0 || i == 1))
	{
//		clog << "END" << endl;
		return cnt.at(i) = i;
	}
…
        /*计算当前位置，保存在变量result中，如果空间不够，还要开辟新空间*/
…
        /* storage ，为返回上一级准备*/
	 if (layer != preview)
	 {
		  storage(curName);
		  download(preName);
	 }
        return result;
}

总体设计完成。接下来就是main函数的编写，以及调试用途的脚手架，详见下文。

运行代码

/* Euler14 Longest Collatz sequence n -> n / 2 , n is even * n -> 3n + 1 n
   is odd * 找出100万下数链最长的n值 * 使用外存 */
#include<cassert>
#include<cstdlib>
#include<ctime>
#include<fstream>
#include<iostream>
#include<sstream>
#include<vector>

using namespace std;

const int LEN = 1000;
const int SIZE = 4 * LEN;
vector < int >cnt(SIZE);
int totalSize = SIZE;			// 记录总空间
stringstream sstream;


void download(const string & name)
{
	cnt.assign(SIZE, 0);
//	clog << "clear" << endl;
	fstream file(name, ios::binary | ios::in);
	if (file)
	{
		assert(file.is_open());
//		clog << "downloading..." << name << endl;
		for (int i = 0; i < SIZE && !file.eof(); i++)
		{
			file.read((char *)&cnt.at(i), sizeof(int));
		}
		file.close();
/*
		for (int i = 0; i < SIZE; i++)
		{
			clog << cnt.at(i) << " ";
		} clog << "
accomplished" << endl;
*/
	}							// endif
}

void storage(const string & name)
{
	fstream file(name, ios::binary | ios::out);
	if (!file)
		throw name;
//	clog << "storaging..." << name << endl;
	for (int i = 0; i < SIZE; i++)
	{
		file.write((char *)&cnt.at(i), sizeof(int));
	}
	file.close();
/*
	for (int i = 0; i < SIZE; i++)
	{
		clog << cnt.at(i) << " ";
	}
	clog << "
accomplished" << endl;
*/
}

int collatz(const int layer, const int i, const int preview)	// preview用于模拟弹栈,表示上一级
{
	assert(i >= 0 && i < SIZE && (layer + 1) * SIZE <= totalSize);
	int n, result;

	/* download */
//	clog << "(" << layer << "," << i << "," << layer * SIZE + i << "," << preview << ")" << endl;
//	system("sleep 1");
	string curName = "./Euler14-data/Euler14-";
	string preName = "./Euler14-data/Euler14-";;
	string suffix;
	sstream << layer << endl;
	sstream >> suffix;
	curName.append(suffix + ".data");
	suffix.clear();
	sstream << preview << endl;
	sstream >> suffix;
	preName.append(suffix + ".data");
	if (layer != preview)
	{
		if (preview >= 0)
		    storage(preName);
		download(curName);
	}


	if (layer == 0 && (i == 0 || i == 1))
	{
//		clog << "END" << endl;
		return cnt.at(i) = i;			/* 边界条件 */
	}
	
	if (!cnt.at(i))
	{
		if (i % 2)
		{
			n = 3 * (layer * SIZE + i) + 1;

			while (n >= totalSize)
			{
//				clog << "懒汉式开辟外存空间..." << endl;
				totalSize += SIZE;
				if (totalSize < 0)
					throw - 2;
			}

		}
		else
		{
			n = (layer * SIZE + i) / 2;
		}
		/* assert cnt[j != i] immutable ] */
		int temp = collatz(n / SIZE, n % SIZE, layer) + 1;
		cnt.at(i) = temp;
//		clog << "Back to " << layer *SIZE + i << endl;
	}							// endif
//	else  clog << "END" << endl;
	result = cnt.at(i);	
/*	
for (int i = 0; i < SIZE; i++)
	{
		clog << cnt.at(i) << " ";
	}
	clog << endl;
*/
	/* storage */
	 if (layer != preview)
	 {
		  storage(curName);
		  download(preName);
	 }
	return result;
}

int main()
{
	unsigned int start = clock();
	try
	{
		cnt.at(0) = collatz(0, 0, -1);
		for (int i = 1; i < LEN; i++)
		{
			if (i / (LEN / 1000) * (LEN / 1000) == i)
			{
				system("busybox clear"); //清理界面
				clog << "[" << i / (LEN / 1000) << "‰]" << endl;
//				clog << "[" << i  << endl;
			}
			cnt.at(i) = collatz(0, i, 0);
			if (cnt.at(i) < 0)
				throw - 1;
		}
//		storage(curName); //保留未完全转换的

		int maxLen = 0, maxI = -1;
		for (int i = 1; i < LEN; i++)
		{
			// cout << cnt.at(i) << " ";
			if (maxLen < cnt.at(i))
			{
				maxLen = cnt.at(i);
				maxI = i;
			}
		}
		cout << "maxLen=" << maxLen << "(" << maxI << ")" << endl;
	}
	catch(string e)
	{
		cerr << e << "打开失败" << endl;
	}
	catch(int e)
	{
		switch (e)
		{
		case -1:
			cerr << "cnt数组存在负数" << endl;
			break;
		case -2:
			cerr << "totalSize整数越界" << endl;
			break;
		}
	}
	unsigned int end = clock();
	cout << "所用时间" << (end - start) / CLOCKS_PER_SEC << "S" << endl;
	return 0;
}

当然，我的运行环境是安卓系统。相关的system调用与电脑略微不同。读者需修改一下才能在计算机上运行。