最近遇到一个题,题目是这样描述的:对于一个数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 << "\naccomplished" << 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 << "\naccomplished" << 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调用与电脑略微不同。读者需修改一下才能在计算机上运行。