如何让学习Wolfram Language更有趣（上）。

Wolfram

NO.1

引言

NO.2

Josephus问题 

我们将编写一个非常简短的程序来计算结果。程序的输入将是40个数字的列表：
{1,2,3,…,40}。Mathematica有一种使用函数Range生成此类列表的简洁方式。在新单元格中输入：Range[40]
然后按住Shift键同时按回车键。尝试在该函数中输入除40以外的其他值。

我们需要的另一个函数是：RotateLeft。如果在新单元格中输入：?RotateLeft，Mathematica将简要说明该函数的功能；如果鼠标单击继续符号“>>”，浏览器将带您前往Mathematica网站，查看带有示例的更详细解释（这对所有Mathematica函数均适用，我们鼓励读者每遇到一个新函数都这样做）。创建新单元格后，尝试计算 RotateLeft[Range[10], 3] 或 RotateLeft[Range[12], 5]。

对于我们的问题，我们只需要 RotateLeft[x]（它等同于 RotateLeft[x, 1]）。思路是“保存”列表开头的那个人员，通过 RotateLeft 将该人员移到列表末尾；在“旋转”之后，我们删除新的列表头部（即被原列表头部刺死的那个人）；这通过函数 Rest 实现。给定一个列表，Rest 函数会删除列表的头部并保留列表的其余部分。在新单元格中尝试执行：Rest[RotateLeft[Range[40]]]。它是否达到了您的预期？我们计划反复使用这种旋转/删除组合，直到只剩下一个人为止。

由于我们可能希望在此处使用除40以外的其他数字，为何不用变量替换 ‘Range[40]’ 呢？也就是说，让我们定义一个函数 Rest[RotateLeft[x]]，打开一个新单元格并输入：f[x_] := Rest[RotateLeft[x]]。然后使用 Shift+Return 将其输入 Mathematica。

首先请注意，定义左侧引入的变量后面必须带有下划线；函数定义中引入的每个变量后都必须跟一个下划线。另请注意使用的是“:=”而不仅仅是“=”。“:=”表示这是一个函数的定义，每当“x_”被替换为可计算的值（例如 Range[40]）时，该函数将重新求值。

计算 f[Range[40]]。注意第一个被淘汰的是谁。接下来计算 f[f[Range[40]]]。我们已经从四十人的圆圈中淘汰了两人；每次应用 f 时，我们都先旋转，然后删除当前队列中的第一个。

现在应该很明显，我们需要将 f 应用于 Range[40] 共三十九次。

当然，我们更希望不必“手动”完成这一操作。（想想有2000人的圆圈！）幸运的是，Mathematica 有一个内置函数可以为我们完成这一任务，即函数 Nest。（在新单元格中输入：?Nest 并回车以获取更多信息）。Nest 让我们可以轻松计算诸如 g[g[g[g[a]]]] 这样的表达式。我们只需告诉 Nest 函数是什么（本例中为“g”）、应用于什么（本例中为“a”）以及应用多少次（本例中为“4”次）。因此，我们可以写 Nest[g, a, 4] 来代替 g[g[g[g[a]]]]。试试看。

我们需要将函数 f 应用于 Range[40] 共三十九次，因此我们要计算：Nest[f, Range[40], 39]。计算该表达式，看看谁幸存下来！

更一般地，我们定义如下函数：survivor[n_] := First[Nest[f, Range[n], n – 1]]。函数“First”用于去掉括号；因此输出不再是“{17}”，而是直接得到“17”。（通常，对列表头部应用 First，会给出列表的第一个元素。）

使用这个 survivor 函数做一些探究：随着圆圈中人数增加，它的表现如何？您能在这里找到规律吗？学习如何使用 Table 函数，以便“批量”获取和显示结果。例如，尝试输入：Table[{n, survivor[n]}, {n, 1, 100}]。

通过使用Mathematica的TableForm函数可以获得更好的格式，例如：TableForm[Table[{n, survivor[n]}, {n, 1, 100}]]。试试看！

您能从这些结果中看出规律吗？如果无法立即识别出规律，对它们“绘图”可能会有所帮助。输入：t = Table[{n, survivor[n]}, {n, 1, 100}];（分号“;”用于抑制输出），然后输入：ListPlot[t]。

练习。尝试对幸存者函数的计算方法提出一个猜想；特别是，您能否手动推算出 survivor[2000] 的值是多少。

NO.3

Josephus

问题的一个更好解法

NO.4

最佳解法 

Josephus问题的一个优雅解法出现在Graham、Knuth和Patashnik合著的《具体数学》（Concrete Mathematics，Addison Wesley，1994）一书中。该解法使用整数的二进制（即以2为底）表示，而非通常的十进制表示。因此，假设 xx 的二进制表示为 bkbk−1bk−2…b2b1b0bkbk−1bk−2…b2b1b0，其中 bk=1bk=1，且每个 bi,0≤i

练习:试用这个新函数，并与之前定义的函数进行比较。

NO.5

3x + 1 问题 

NO.6

素数

们将编写一个非常简短的程序来计算结果。程序的输入将是40个数字的列表：

{1,2,3,…,40}。Mathematica有一种使用函数Range生成此类列表的简洁方式。在新单元格中输入：Range[40]

然后按住Shift键同时按回车键。尝试在该函数中输入除40以外的其他值。

练习。 证明当 0≤n≤390≤n≤39 时，这个二次式能生成40个素数。对于更大的 nn 值呢？提示：你可以通过 (#2+#+41)&/@Range[0,39](#2+#+41)&/@Range[0,39] 轻松生成待测试素性的列表，其中 f/@wf/@w 表示将函数 ff 应用于列表 ww 的每个元素后得到的列表。另一种方法是使用 Table：Table[i^2 + i + 41, {i, 0, 39}] 来生成该列表。

接下来我们研究正整数中素数的“密度”。我们考察区间 [1,…,n][1,…,n] 中素数的个数与 nn 的比值。Mathematica 有一个内置函数 PrimePi 用于计算该区间内素数的个数，我们将使用它。例如，要找出小于 1000 的素数个数，

或者，我们可以使用 Length[sieve[n])。

我们不打算计算该函数的值，而是使用 Plot 命令对其绘图。

这个图让19世纪的数学家联想到了1/Log[x].的图形。

事实上，可以证明当 x→∞x→∞ 时，prdensity[x] 与 1log⁡[x]log[x]1 的比值趋近于1，这一重要结果被称为素数定理。

练习 :对越来越大的 x 值计算该比值。

孪生素数

孪生素数是指相差为2的一对素数；例如，(5,7)是孪生素数，(17,19)也是。我们希望找出小于某个数的所有孪生素数。首先使用上面定义的函数 primelist 列出区间{n,m} 内的所有素数。

然后使用 Select 结合 MemberQ，从 primelist 中选取满足 #+2 也属于 primelist 的元素 #。

接下来，我们展示如何从这些列表重新构造孪生素数对。

我们希望列出这两个列表中对应元素组成的数对。使用 Mathematica 命令 MapThread。给定两个列表 {a,b,c} 和 {d,e,f} 以及一个函数 h，命令 MapThread[h, {{a,b,c},{d,e,f}}] 将 h 应用于这两个列表的对应元素。

这里我们需要函数h[x,y]={x,y}。

数学中最著名的未解问题之一是孪生素数猜想，它断言存在无穷多对孪生素数。

练习：找一些大的孪生素数对。

练习：将素数之间的“间隔”2替换为其他偶数，并在大区间中找到一些这样的数对。这些数对的“密度”与孪生素数的密度相比似乎如何？

梅森素数与完全数

形如 2n−12n−1 的素数称为梅森素数，以17世纪研究它们的法国神父马林·梅森的名字命名。

练习: 找出所有小于1000的梅森素数。提示：首先生成小于 210210（=1024=1024）的形如 2n−12n−1 的数；然后从该列表中选出素数。

如果一个正整数等于其所有正真因子之和，则称它为完全数。注意，数本身不是它自己的真因子。因此6是完全数，因为 6=3+2+16=3+2+1。同样，28也是完全数，因为 28=14+7+4+2+128=14+7+4+2+1。要找出所有小于1000的完全数，我们可以使用 Mathematica 的测试函数 PerfectNumberQ。此外，Mathematica 函数 PerfectNumber 会给出第 n 个完全数。

如果 pp 是梅森素数，则 p(p+1)/2p(p+1)/2 是完全数。这是欧几里得证明的！事实上，欧拉证明了所有偶完全数都具有这种形式。我们可以使用 Mathematica 的 Solve 函数轻松验证这一结果。例如，以下代码验证了欧拉关于 8128 的结果。

练习:找一些更大的完全数并验证欧拉的结果。

目前尚不清楚是否存在无穷多个完全数，也不清楚是否存在奇完全数。

然而，欧几里得证明了，只要 (2ρ−1)(2ρ−1) 是（梅森）素数，则 2ρ−1(2ρ−1)2ρ−1(2ρ−1) 总是一个偶完全数。

欧拉证明了所有偶完全数都具有这种形式。

练习：对前10个完全数验证欧拉的结果。

索菲·热尔曼素数

若素数p 满足2p+1 也是素数，则称p 为索菲·热尔曼素数。

练习： 找出所有小于 10,000 的索菲·热尔曼素数。

整数分解为素数

Mathematica 可以使用命令 FactorInteger 轻松地将整数分解为其素数因子。结果是一个数对列表，其中第一个元素是素数，第二个元素是该素数在分解中应取的幂次。

我们也可以使用 FactorInteger 来分解一个数字列表。如果你查询 FactorInteger（?FactorInteger），你会在“范围”下看到，它会逐元素地在线程化方式下处理列表。这允许进行以下计算。然而，如果你需要分解非常大的整数，请准备好等待很长时间！

哥德巴赫猜想

数学中一个非常著名的未证明猜想是哥德巴赫猜想，它断言每个大于4的偶数都可以写成两个奇素数之和。我们可以使用 FindInstance 来研究它。因此，要证明这对582成立，我们可以计算以下内容。

练习：使用 FindInstance 证明哥德巴赫猜想对于不超过600的偶数成立。

对于偶数n≥670，FindInstance 无法返回结果。

不过，还有另一种方法可以验证哥德巴赫猜想对不超过某个界限的偶数成立。思路是从一个初始素数段开始，如下所示。

然后，我们使用 Tuples[ , +2] 从该列表生成所有二元组列表。

接着，我们使用 Plus 对这些二元组求和，并将其线程化（threading）应用于这些列表。我们使用的命令是

Apply[Plus, Tuples[pr[50], 2], {1}] 或其简写形式 Plus@@@Tuples[pr[50], 2]。（查阅 Apply。）最后，我们在其前面加上 Union 以消除重复项，从而形成最终列表，该列表包含能表示为两个素数之和的所有数。

注意，在列表末尾，一些偶数缺失；推测它们会在后续列表中出现。更一般地，我们在下面定义函数 goldbach。

练习： 使用函数 goldbach 研究哥德巴赫猜想。特别地，证明 670 满足该猜想。

不必检查特定数字的列表，而是使用 MemberQ[List, a] 来表明 a 是否属于 List。

你能验证哥德巴赫猜想的最大偶数是多少？

欧拉乘积公式

1737年，L·欧拉证明了一个涉及无穷和、无穷乘积以及素数的非凡公式。

首先我们介绍Mathematica中的求和与求积函数。Sum 和 Product 是每个Mathematica用户都需要知道的函数。要求前100个正整数之和，可以输入：

如果我们想要前 nn 个正整数之和，可以输入：

因此，得到了这个和式的著名公式。如果我们使用面板中的“基本数学助手”，可以按照如下方式使用标准的 Σ- 符号：

练习。 使用 Mathematica 找出：

的等式。

类似地，我们可以使用 Product 或 ∏∏。

在 Sum 或 Product 中，上界可以是 Infinity（或 ∞）。因此，

最后一个结果是下面定义的黎曼ζ函数的一个特例。

最后一个结果是下面定义的黎曼ζ函数的一个特例。

我们也可以定义涉及无穷乘积的函数，例如，其中 Prime[i] 是第 i 个素数。

练习： 证明上面定义的两个函数 ζ,ηζ,η 在 s=2,…,100s=2,…,100 上相等。提示：使用 Table 和 == 来验证这一点。
这一般性结果，即这两个函数对所有 s>1s>1 相等，由 L. Euler 于1737年给出，被称为黎曼ζ函数的欧拉乘积公式。