汉诺塔的五种解法

· · 算法·理论

0x00 前言

汉诺塔是一个很古老的问题了，本文将介绍求解该问题的五种实现方式。

汉诺塔有五种解法，你知道么？ ——鲁迅

0x01 问题介绍

汉诺塔问题简单来说就是给三根柱子，一开始其中一根上叠放着多个圆盘，圆盘由上到下直径依次增大，另外两根柱子为空。按一定的规则在柱子间移动圆盘，求将所有圆盘移动到另一根柱子上的最小步数。

移动规则：

0x02 解法

如果记圆盘数量为 n，对应的最小移动步数为 f_n，则显然有

\begin{aligned} f_1&=1\\ f_2&=3\\ f_3&=7\\ &\cdots \end{aligned}

（自己想一想就明白了）

但是如果 n 比较大，我们就需要 f_n 的通项公式了。

其实从上面的列表就可以猜出 f_n=2^n-1。

事实上这是对的，证明如下：

首先求递推式：

标记三个柱子分别为 A,B,C，圆盘一开始在 A 上，最终要挪到 C 上。假设我们已经求出了把 n-1 个圆盘从一个柱子挪到另一个柱子上上需要的 f_{n-1} 步，则将 n 个圆盘从 A 挪到 C 上，相当于以下三部分的工作：

  1. 用 f_{n-1} 步将 A 最上面的 n-1 个圆盘挪到 B 上
  2. 用 1 步将 A 上剩的最大的圆盘挪到 C 上
  3. 用 f_{n-1} 步将 B 上的 n-1 个圆盘挪到 C 上

由于挪较小的 n-1 个圆盘时最大的圆盘不会干扰，以上三步是成立的。

总步数 f_n=2f_{n-1}+1。

利用高中的数列知识，可以求出 f_n 的通项：

f_n=2f_{n-1}+1 f_n+1=2(f_{n-1}+1) \because f_1=1 \therefore f_n+1=2^n f_n=2^n-1 \mathrm{Q.E.D.}

然后用高精度你就能做 P1760 了。

0x03 编程解决

以下默认是模意义下的解决办法，且 #define int unsigned long long，const int mod = 1e9+7。

代码均在 Xcode 内，C++14(GCC9) 下编写并运行。

如果仅使用递推式，我们有两种解法。

1. 递归

前置知识：C++ 语言基础。

int hanoi1(int n) {
    if (n == 1) return 1;
    return (2 * hanoi1(n - 1) + 1) % mod;
}

时空均 O(n)，容易爆栈 \texttt{MLE}。

2. 递推

前置知识：C++ 语言基础。

int hanoi2(int n) {
    int x = 1;
    while (n--) x = (x * 2 + 1) % mod;
    return x;
}

时间 O(n)，空间 O(1)，1\text s 内处理 10^8 的 n 没有问题。

但是如果 n 更大，我们就不得不在利用通项 2^n-1 的同时，使用一些神奇的算法了。

3.快速幂

前置知识：幂的性质。

以下内容是为没学过快速幂的萌新准备的，大佬请跳过。

快速幂就是用来在 O(\log n) 的时间内处理某一个数自乘 n 次的算法。

例如，给出 a，要计算 a^{11}。

最简单的就是循环暴力乘，但是这样是 O(n) 的，很慢。

但是如果我们把 11 写成二进制：11=(1011)_2

即：11=2^3+2^1+2^0

再利用幂的性质：a^{x+y}=a^xa^y

就可以得出 a^{11}=a^{2^3}a^{2^1}a^{2^0}。

乍一看好像变复杂了，但是我们发现，a^{2^i} 可以通过 i 次 a=a*a 的循环得出，而最大可能的 i 是指数的二进制的最高位，是 O(\log n) 的。

又注意到一个 i 出现在指数上，当且仅当 n（这里是 11）的二进制的第 i 位（从 0 开始）是 1。

所以可以写出如下代码：

int quickpow(int a, int n) {//a^n%mod
    int ans = 1;
    while (n) {
        if (n & 1) ans = ans * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return ans;
}

用来过模板题。

而对于汉诺塔，底数 a=2，最后还要减一，于是代码变成：

int hanoi3(int n) {
    int ans = 1, a = 2;
    for (int i = 1; i <= n; i <<= 1) {
        if (n & 1) ans = ans * a % mod;
        a = a * a % mod;
    }
    return ans - 1;
}

时间 O(\log n)，空间 O(1)，对于 [0,2^{64}-1] 范围内的数据都能给出答案。

4. ???

前置知识：位运算。

我们要求的是 2^n-1。

回想一下，位运算中 1<<n 表示的就是 2^n。

但是 n 太大了，会爆炸。

所以将 n 写成 aL+b，其中 L 视情况定为不会使 1<<L 爆炸的数，而 b<L，这样就可以 O(1) 求 x=2^L,y=2^b，再快速幂求 z=x^a，最后得到答案为 yz-1。

代码：

int hanoi4(int n) {
    const int L = 63;
    int a = n / L, b = n - a * L;
    int x = 1, y = 1, z = 1;
    x = (x << L) % mod;
    y = (y << b) % mod;
    while (a) {
        if (a & 1) z = z * x % mod;
        x = x * x % mod;
        a >>= 1;
    }
    return (y * z - 1) % mod;
}

时间 O(\log n)，空间 O(1)，但较方法 3 有因数而异的常数优化。

当然可以对该方法进行各种奇奇怪怪的修改，不幸的是它们的常数都超过了原始算法。

以上都是对单次询问的解法。如果有多次询问，则又有更优的做法。

5.光速幂

前置知识：与约数有关的知识、模算术。

其实方法 4 就借用了光速幂的思想。

对于一般的底数：

假设模数为 $p$，预处理 $a^1,a^2,\cdots,a^L,a^{2L},\cdots,a^{\lceil\frac p L\rceil\times L}$ 模 $p$ 的结果，则 $a^{kL}$ 和 $a^b$ 都在其中，就可以 $O(1)$ 计算 $a^n$ 了。 但是为什么处理到 $a^{\lceil\frac p L\rceil\times L}$ 就可以了呢？ 这里要用到一个叫欧拉定理的东西。 这个定理先定义了一个欧拉函数 $\varphi(n)$，表示小于正整数 $n$ 的正整数中与 $n$ 互质（最大公约数为 $1$）的个数。例如：小于 $10$ 的正整数中有 $1,3,7,9$ 四个与 $10$ 互质，所以 $\varphi(10)=4$。 欧拉定理：若正整数 $a,n$ 互质，则 $a^{\varphi(n)}\equiv1\pmod n$。 容易看出，在 $n$ 是质数时，$\varphi(n)=n-1$。 （题外话：显然费马小定理就是欧拉定理在 $n$ 是质数时的形式。） 欧拉定理的证明就不讲了，有兴趣的可以自己去找。 在此 $\texttt{\%\%\%Euler}$！ 然后我们要用到的，是扩展欧拉定理： $a^b\equiv a^{b\bmod\varphi(p)}\pmod p

其实很简单，就是 a^b\equiv a^{\varphi(p)\lfloor\frac b{\varphi(p)}\rfloor+b\bmod\varphi(p)}\equiv 1\cdot a^{b\bmod\varphi(p)}\pmod p。

根据定义，有 \varphi(p)\le p-1，所以指数 b\bmod\varphi(p)\le\varphi(p)-1\le p-2。

对于之前的问题 a^n，

n=kL+b\le p-2 kL\le p-2-b\le p-2 所以处理到 $\lceil\frac p L\rceil$ 就可以了。 那么 $L$ 取多少呢？ 来算算预处理时间复杂度：$O(L+\frac p L)=O(\max(L,\frac p L))

最小值发生 L=\sqrt p 处，所以取 L=\sqrt p 最优。

预处理代码：

int f[32000], g[32000];//分别存a^1~L,a^L~p/L*L
int a = 2;//底数
int L = sqrt(p);//sqrt函数需要#include <cmath>
f[0] = g[0] = 1;
for (int i = 1; i <= L; i++)
    f[i] = f[i-1] * a % p;
for (int i = 1; i <= L; i++)
    g[i] = g[i-1] * f[L] % p;

然后求模数的欧拉函数，本来应该讲讲的，但是我懒模数 =10^9+7 是质数，所以欧拉函数就是模数减一。

最后，2^n 就是 f[n%L] * g[n/L] % mod。

所以

int f[32000], g[32000];
int L;
void init() {
    L = sqrt(mod);
    f[0] = g[0] = 1;
    for (int i = 1; i <= L; i++)
        f[i] = f[i-1] * 2 % mod;
    for (int i = 1; i <= L; i++)
        g[i] = g[i-1] * f[L] % mod;
}
int hanoi5(int n) {
    n %= mod - 1;
    return (f[n%L] * g[n/L] - 1) % mod;
}

用的时候 O(\sqrt p) init() 一遍，然后就可以疯狂 O(1) 计算了。

空间复杂度 O(\sqrt p)。

模板题

0x04 各解法运行时间比较

使用代码：

#include <cstdio>
#include <ctime>
/*...hanoi1~5...*/
void test(int n, int t) {
    clock_t t0, t1, t2, t3, t4, t5;
    t0 = clock();
    for (int i = 0; i < t; i++) hanoi1(n);
    t1 = clock();
    for (int i = 0; i < t; i++) hanoi2(n);
    t2 = clock();
    for (int i = 0; i < t; i++) hanoi3(n);
    t3 = clock();
    for (int i = 0; i < t; i++) hanoi4(n);
    t4 = clock();
    init();
    for (int i = 0; i < t; i++) hanoi5(n);
    t5 = clock();
    printf("1: %lu\n", t1 - t0);
    printf("2: %lu\n", t2 - t1);
    printf("3: %lu\n", t3 - t2);
    printf("4: %lu\n", t4 - t3);
    printf("5: %lu\n", t5 - t4);
}

多次结果取平均后整理如下表（注意 \#5 最后一列是 \times10^3）：

(n,t)~ ~(10^5,1000)~ ~(10^8,10)~ ~(100,5\times10^7)~ ~(10^{18},10^6)~
\#1 7.18\times10^5 MLE TLE MLE
\#2 4.34\times10^5 4.40\times10^5 TLE TLE
\#3 51 6 6.86\times10^5 2.88\times10^5
\#4 32 1 2.39\times10^5 2.70\times10^5
\#5 440 440 1.32\times10^5 3.13\times10^3

0x05 结语

是不是感觉费这么大的劲，就做了道红题？

但是，探索解决同一问题的不同方法总是好的。这可以帮助我们巩固已学知识，并在知识之间建立联系，加深对知识的理解，从而引发思考。

总之，受益匪浅。

希望各位在今后的学习中，也要勤于思考，敢于创新。

0x06 参考文献