因为偶然的机会听说了An Applied Mathematician’s Apology一书，我对此很感兴趣就去读了一下，收获不少。这本书的作者是Lloyd N. Trefethen教授，一位在数值计算领域相当著名的老教授，我的数值线性代数课程所用的教材就是他写的。书中的许多观点都让我能感到共鸣，但也有些地方让人感觉到“时代变了”，与当今时代科研的关注点并不尽相同。
首先，最重要的是：计算数学是一门实验科学。这是Trefethen教授在书中多次强调的观点，也让我相当认同。无论是教授所述的经历还是我自己的经验，科研发现总是来自于对程序运行结果、特别是可视化结果的观察，而不是来自数学推导。和计算机科学一样，计算机是我们不可或缺的实验工具。
关于计算的背后是否需要严谨的数学理论，我的理解一直在改变：对于OIer时期的我来说答案肯定是“不需要”，而从高三到大二我的观点逐渐转向了“需要”，而现在我的回答会是“需要，但前提是理论假设与实际数据足够契合”。在数值分析中，我们经常关注算法的收敛性和稳定性，这当然是有意义的；费了九牛二虎之力编写了计算程序却发现它对于实际数据压根就不收敛，或者一加噪声结果质量就急剧恶化，结局就是白干。但是这种分析并不代表一切，比如Gauss消元就算有pivoting，误差理论上也是指数增长的，但实际上这是一个很有效的算法。这是为什么呢？因为大自然并不是一个adversary，它不会存心来hack你的数值计算程序，我们关心的是对于实际数据表现如何，而不是理论最坏。对此，Trefethen教授在该书中有一句评价让我相当认同：
If rounding errors vanished, 90% of numerical analysis would remain.
人们对于数值分析总是太过悲观了，认为它是一个丑陋的应用学科，总是在应对数据和运算的不完美。而Trefethen教授认为，数值分析应该被定义成“连续数学中的算法研究”，正如计算机科学是“离散数学中的算法研究”一样。
在纯数学的analysis中，我们关心的总是两件事：
- 空间的性质：有限维还是无穷维？在$\mathbb{R}$上还是在$\mathbb{C}$上？配备什么norm？完备吗？可分吗？$\dots$
- 函数的性质：有界，连续，可导，解析，$\dots$
这两件事确实在计算中是至关重要的，也经常是对数值算法的理论分析中经常用到的。但实际上，这些数学的分析与实际数据并不契合。以求解方程$f(x)=0$为例，如果$f$足够光滑，直接Newton迭代就能在$O(\log \log \epsilon^{-1})$的时间之内求出解了；而如果只把$f$看成没有任何特殊性质的函数，显然除了逐个数据点枚举之外没有任何有效的方法。而对于现实中遇到的函数来说，这两种假设都太过极端了。因此，我认为应该把这两点做如下的替换：
- 空间的性质：低维还是高维！对于underlying field是什么我并不关心，无论是实数还是复数，还是实际运算中的浮点数，抑或是图片中的RGB三元组，对于计算来讲这都不重要。无限维？在计算中这是和实数一样荒唐的假设。但低维和高维之间是真的有本质区别的。
- 函数的性质
- 对于低维的数据，Trefethen教授在书中的总结相当精辟：分段光滑。不是悲观的$L^\infty$，也不是乐观的$C^\infty$，忘掉这些没用的假设吧。从全局来看，实际的数据没有任何好的性质；而从局部来看，可以认为是光滑的。一个很直观且常见的例子就是此时此刻你眼前的这幅二维图像，它就可以被视为一个贴上了texture的分段光滑函数。视频或音频数据也是如此，数值计算中的PDE求解也是如此，FMM也利用了这一点。
- 而高维的数据就完全不同了，迫于时代的局限Trefethen教授在书中完全没有谈论这一部分，而这在今日的世界中越来越重要。如果我们孤立地看一幅图片，它是一个二维的数据；而如果我们使用data-driven的思维来看，考虑所有的人脸图片构成的数据集，这就是一个高维的数据了。对于高维的数据，函数就是没有任何性质的，任何精妙的算法都会被curse of dimensionality碾成渣，我们只能把所有映射一视同仁看成机器学习的黑箱，所能依赖的只有monte carlo和堆算力。
这个话题就有些扯远了，下面我回到本文的主题上来：Ten Digit Problems。这是这本书中最让我感兴趣的部分，Trefethen教授给他的博士生提出了一些有趣的数值计算小问题，并要求把答案计算到10位有效数字。为什么是10位呢？虽然大部分实际应用中我们并不需要如此高的精度，但很多时候我们想要解决的问题（如求矩阵的特征值）是作为一个模块被实际应用的程序调用的，这种时候往往就会有更高的精度要求了，而10位精度往往是使用64位浮点类型能达到的合理精度要求。虽然要求是10位，其实Trefethen教授希望把答案计算到尽量多位的数字，比如说10000位。有一篇文章更全面地介绍这些Ten Digit Problems的前因后果。而在我看来，这种高精度计算就和数值计算关系不大了，而是属于纯数学和CS的领域：你需要先使用高端的数学技巧写出答案的封闭形式、牛顿迭代式或者超几何级数表达，然后使用Binary splitting在超级计算机上算到几亿位小数。在下文中我只会关注如何求解An Applied Mathematician’s Apology一书中列出的那些问题，达到10位有效数字，并且和AI对比（使用的AI是GPT-5.2 Thinking）。两个微分方程的题我就跳过了，因为我对此不太了解也不感兴趣。
第一题：粒子到达底行的概率
A particle starts at the top vertex of a triangular array with 30 points on each side, and then takes 60 random steps. What is the probability that it ends up in the bottom row?
额$\dots$书中Trefethen教授提到他虽然和鼎鼎大名的Knuth教授属于同一个部门（Stanford CS），但对于Knuth教授的研究完全不了解，看来所言非虚啊。总之，对于任何一个有点CS基础的人来说，这个题都是极为简单的，和其它题完全不属于同一个级别的。但毕竟要考虑到时代的变迁，在1994年数字三角形还是IOI题呢。
解法就不多说了。对于这么小的数据规模，直接暴力递推即可。我秒了，意料之中AI也秒了。C++代码，答案是
\[\boxed{9.512343502\times 10^{-6}.}\]第二题：sin(n)/log(n)的求和
What is $\sum_{n=2}^{\infty} \sin(n)/\log(n)$?
这道题对我来讲就困难很多了，属于一个我不太了解的领域，但是思考了一会之后仍然成功解决了。它的主要难度在于数列$\{\frac{\sin(n)}{\log(n)}\}$的收敛实在是太慢了，直接暴力求和当然是行不通的。如果要求的不是求和而是积分的话，我们会如何加速呢？一种想法是换元，但是换元对于求和来讲基本是个不存在的操作；另一种想法就是分部积分：给$\sin x$积分而给$\frac{1}{\log x}$求导，由于$\sin x$在积分过程中仍然是有界的，而
\[\frac{d}{dx^k}(\frac{1}{\log x})=O(x^{-k}\log^{-2} x),\]分部积分几次之后收敛就会很快了。那么求和也能“分部求和”吗？确实是可以的，而且是我在最基本的数学分析课程中就学过的：Abel transformation. 在具体讲解之前我们需要先约定一套离散的微积分记号：记$\Delta$为差分算子，
\[\Delta[f](x)=f(x+1)-f(x).\]如果$\lim_{n\to\infty}f_ng_n=0$，就有如下的公式：
\[\sum_{k=m}^{\infty}f_k\Delta g_k=-f_{m}g_{m}-\sum_{k=m}^{\infty}g_{k+1}\Delta f_k.\]该公式很容易证明，首先证明有限求和的版本，然后取极限即得证。在实际操作中，我们对$e^{in}$做等比数列求和，而对$\frac{1}{\log n}$做差分，在分部求和5次之后暴力求和即可。C++代码，答案是
\[\boxed{0.6839137864.}\]AI轻松解决了这个问题。一开始它试图调用mpmath.nsum直接逃课，而我禁止了逃课之后它也很容易就解决了，通过Abel summation，即给$\sum a_n e^{-nx}$求和再取$x\to 0^{+}$的极限。这其实是我一开始的想法，但是我害怕外推是数值不稳定的，就没有下手去做，而是换了一个想法。AI使用Neville 插值算法实现的外推，并且发现用8个点的外推结果就足够好了。
第三题：正四面体包装问题
Three regular tetrahedra each have volume 1. What’s the volume of the smallest sphere you can fit them inside?
这个问题很类似Circle packing，本质上是极为困难的，但是这里的数据范围足够小。一个naive的猜测是三个正四面体共用一条棱的时候取到最小值，此时答案为$3\sqrt{6}\pi\approx 23.0858969429136$；虽然这个配置很对称，看上去填充率也不错，但这并不一定代表它是全局的最小值；实际上还是需要先寻找一个最优的配置，再针对这一配置精确求解。这也是我们之后会经常用到的一个思想：首先找到精度很低的近似解，再利用这个近似解的性质去高精度求值。
如何低精度地寻找一个最优的配置呢？实际上我的第一篇博客就是讨论Circle packing问题的，我早已了解模拟退火对此很有效。当然有更好的优化算法，但是对于这么小的规模，模拟退火就很够用了。说起来简单做起来难，虽然我在使用有限种颜色逼近图像一文里已经实现过了3维的随机增量最小球覆盖，但计算正四面体的顶点并判断是否相交是个工作量很大的事情。我虽然会写，但是对于这种纯粹的体力活完全懒得去写，就让ChatGPT给我生成了模拟退火的代码。这么复杂的任务AI是没法一次写对的，但是我给了它我的最小球覆盖代码之后它就成功完成了，代码见p3_sa_chatgpt.cpp。
模拟退火的结果显示我之前的naive猜测确实错了，它跑出来的解体积约为$22.141631408067$。我又让AI写了一个代码绘制这个解的3维图片，大概是这样的：
可以看见这个解的特点是两个正四面体有一面完全贴合在一起构成双三角锥，而剩下一个正四面体放置在双三角锥的其中一个面上。换一个视角之后其性质就更加明显了：
我们建立一个三维的直角坐标系，而这张图表示的是在XY平面上的投影。可以看到图中过$C,E,X$这三点的圆即为最小外接球的大圆，我们需要找到让外接圆最小的$X$点。为方便起见，记棱长$a=\sqrt[3]{6\sqrt{2}}$，则$A,B$两点的坐标为$(0,0,\pm a/2)$。由于$ABC$是正三角形，$C$点的坐标是$(-\frac{\sqrt{3}}{2} a,0,0)$。正三角形的高是$\frac{\sqrt{6}}{3} a$，而$D$在底面上的投影是中心，在中线的$2/3$处，因此$D(-\frac{\sqrt 3}{6} a,\frac{\sqrt 6}{3} a, 0)$。点$E$的坐标也是确定的，但其实无需算出。点$X$的$y$坐标是$-\frac{\sqrt 6}{3} a$，而显然$z$坐标取0是最优的，唯一的未知值是其$x$坐标。这就是一个标准的初中平面几何问题：
给定平面上的$C,E$两点，如何找到在直线$y=-\frac{\sqrt 6}{3} a$上让三角形$CEX$外接圆半径最小的点$X$？
这算是初中时的我最最喜欢的问题之一了，我当初是自己想出来的，后来也在Matrix67的博客上看到过。由于正弦定理，这等价于最大化$\angle CXE$。而考虑拉格朗日乘数法的几何意义，其必要条件是“$\angle CXE$为一定值的$X$构成的曲线”和直线相切。换句话说，圆$CXE$和直线$y=-\frac{\sqrt 6}{3} a$相切。当然这只是必要条件，会解出来两个解，需要把直线$CE$左侧的那个舍掉。我们设圆心坐标为$(x,y)$，则$(x,y)$位于$CE$的垂直平分线上。我们不需要计算$E$的坐标，只需要$CE$中点的坐标就能算出这个垂直平分线了；而$CE$的中点就是面$ABD$的中点$(-\frac{\sqrt{3}}{18} a, \frac{\sqrt{6}}{9} a, 0)$，于是垂直平分线的方程为
\[y=-2\sqrt{2} (x+\frac{\sqrt{3}}{18} a)+\frac{\sqrt{6}}{9} a.\]令$(x,y)$到$C$的距离和到直线$y=-\frac{\sqrt 6}{3} a$的距离相等：
\[\left(x+\frac{\sqrt{3}}{2}a\right)^2+\left(-2\sqrt{2}\left(x+\frac{\sqrt{3}}{18}a\right)+\frac{\sqrt{6}}{9}a\right)^2 = \left(-2\sqrt{2}\left(x+\frac{\sqrt{3}}{18}a\right)+\frac{\sqrt{6}}{9}a+\frac{\sqrt{6}}{3}a\right)^2,\]这是一个$x$的二次方程，解得$x=(-11\sqrt{3}+6\sqrt{10})a/6$。于是外接圆的半径是$4(\sqrt{6}-\sqrt{5}) a$，代入$a$的值就得到体积为
\[512\sqrt{2}(\sqrt{6}-\sqrt{5})^3\pi\approx \boxed{22.1131674629738.}\]顺便扯一下，由于这只是一个二次方程，点$X$是可以尺规作图得出的。这个问题的另一个描述是计算以$C$为焦点、$y=-\frac{\sqrt 6}{3} a$为准线的抛物线和$CE$的垂直平分线的交点，但八年级的我并不知道该如何尺规作图求抛物线和直线的交点（除了暴力模拟解方程），于是向我的同班同学crz请教，他研究了几节课之后还真的找到了一个优美的利用相似三角形的做法，我后来就是使用这个做法通关的Euclidea上的对应关卡。
那么问题的原作者Trefethen教授是如何解决的呢？见MathOverflow。Trefethen教授只是通过摆弄几个正四面体来观察最优情况，而并没有利用优化算法来求出最优配置；而且他也没有解析求解最优情况对应的几何问题，而是数值求解，再使用反符号计算器（如ries，也是我中学时期最感兴趣的东西之一）猜出的解析解。我认为我还是要技高一筹的。
AI的表现就很让我失望了。首先AI试图玩文字游戏逃避问题，直接把三个正四面体重叠在一起；我加上不允许重叠的限制之后，它也完全不知道该如何求解。它编写的搜索效率很差，并不能找到足够好的解；然后AI和最开始的我一样假设该问题具有对称性，并给出了$3\sqrt{6}\pi$这个错误答案。我告诉它不对之后，它仍然在保持这些对称性的前提下搜索；我告诉它并不具有任何的对称性，AI就完全不知道该如何搜索了，它试图爬山，但是效率很差爬不到解；也试过模拟退火，但这个问题还是太复杂了，它写不对搜索程序。于是我只能手动告诉它这个近似解，让它针对近似解去优化；AI并不会观察最优解的几何性质，而直接数值优化也做不对，我怎么提点它都没有用。总之在这个问题上我是完胜，看起来复杂但其实整个过程都比较顺利，我的用时甚至在这套题里算少的。
第四题：震荡函数的积分
What is $\int_{0}^1 \sin^2(\tan(\tan(\pi x))) dx$?
这个数值积分的主要困难在于$\sin^2(\tan(\tan(\pi x)))$是一个震荡极为剧烈的函数，想要直接积分极为困难；所幸对于积分我们是有很多的处理办法的，一层一层地换元就能消除全部的震荡。
首先是最内层的$\tan(\pi x)$，在$x$从0走到1的过程中相当于是从$-\infty$走到了$\infty$，因此换元$u=\tan(\pi x),\ dx=\frac{du}{\pi(u^2+1)}$即变为了
\[\frac{1}{\pi}\int_{-\infty}^{\infty} \frac{\sin^2(\tan u)}{u^2+1} du.\]由于$\tan u$的周期是$\pi$，我们把所有的$\mod \pi$相同的$u$合并到一起求和：
\[\frac{1}{\pi}\int_{-\pi/2}^{\pi/2} \sin^2(\tan u )(\sum_{k\in \mathbb{Z}} \frac{1}{1+(u+\pi k)^2}) du.\]括号内部的求和$\sum_{k\in \mathbb{Z}} \frac{1}{1+(u+\pi k)^2}$其实是有封闭形式的。虽然这很不数值计算，但我碰巧知道如何完成这个求和：Poisson summation。具体来说，我们首先考虑$u=0$的情况，即没有偏移量的求和$\sum_{k\in \mathbb{Z}} \frac{1}{1+\pi^2k^2}$。这个形式一看就很像傅里叶变换的结果，让我们尝试一下：$e^{-|x|}$的傅里叶变换为$\frac{2}{1+4\pi^2\xi^2}$，因此这个求和中的$\frac{1}{1+\pi^2k^2}$就是$e^{-2|x|}$的傅里叶变换，用Poisson求和公式就得到
\[\sum_{k\in \mathbb{Z}}\frac{1}{1+\pi^2k^2}=\sum_{k\in \mathbb{Z}}e^{-2|k|}=\frac{1}{\tanh{1}}.\]下面推广到一般情形：$u+\pi k=\pi(u/\pi+k)$，因此相当于原函数平移了$u/\pi$，傅里叶变换就需要乘上$e^{2\pi i\xi u/\pi}=e^{2 i u \xi}$，即
\[\sum_{k\in \mathbb{Z}}\frac{1}{1+(u+\pi k)^2}=\sum_{k\in \mathbb{Z}}e^{-2|k|+2 i u k}=\frac{\sinh(2)}{\cosh(2)-\cos(2u)}.\]我们记这个封闭形式结果为$F(u)$。此外，还可以使用Proofs from THE BOOK中Cotangent and the Herglotz trick一节的结果来计算这个求和；这是我在该书中最喜欢的一节，强烈推荐大家阅读。概括地说，我们可以把$\frac{1}{1+(u+\pi k)^2}$拆分为两个$\frac{1}{ak+b}$之和，其求和为余切函数。但这一方法更麻烦一些，而且收敛性上也有点细节问题，我觉得还是泊松求和更好。
下面该去除第二层的$\tan$函数了。由于$u\in (-\pi/2,\pi/2)$，换元$v=\tan u$之后$v$的范围就是$(-\infty,\infty)$，且$du=\frac{dv}{1+v^2}$。原式变为
\[\frac{1}{\pi}\int_{-\infty}^{\infty} \sin^2 v \frac{F(\arctan v)}{1+v^2} dv.\]仍然是把所有的$\mod \pi$相同的$v$合并到一起求和：
\[\frac{1}{\pi}\int_0^{\pi} \sin^2 v (\sum_{k\in \mathbb{Z}} \frac{F(\arctan(v+\pi k)}{1+(v+\pi k)^2}) dv.\]到此为止我们就消除了全部的震荡，可以数值积分了。这个无限求和太过复杂，因此我没有尝试算出封闭形式，而是数值求和。由于$k\to\infty$时$\frac{F(\arctan(v+\pi k))}{1+(v+\pi k)^2}=O(k^{-2})$，为了算出10位有效数字需要高达$10^{10}$项，显然是不现实的（以这个问题被提出的年代的算力来衡量）。但使用一个简单的拆项，收敛就能大幅提速了：
\[\frac{F(\arctan(v+\pi k)}{1+(v+\pi k)^2}=\frac{F(\arctan(v+\pi k))-F(\pi/2)}{1+(v+\pi k)^2}+\frac{F(\pi/2)}{1+(v+\pi k)^2}.\]后一项的求和就是$F(\pi/2)F(v)$；而泰勒展开给出$F(\arctan(v+\pi k))-F(\pi/2)=O(k^{-2})$，因此前一项无穷求和的单项是$O(k^{-4})$的，算出10位有效数字就只需要$10^{10/3}$项了，成为了完全可以承受的计算。当然前面$F(u)$也可以使用类似的方法快速求和，因此这道题是可以不用任何高端数学技巧只用数值计算完成的。
我的C++代码，其中数值积分使用的是复合五点Gauss-Legendre求积，是直接复制的我数值分析课程作业的代码。答案为
\[\boxed{0.3909921622.}\]这个问题我被AI完爆了，对于AI来说纯数学的复杂求和与复杂积分根本不是事，它直接给出了封闭形式的答案：
\[\boxed{\frac{1-\exp(-2\tanh(1))}{2}.}\]具体来讲，我们不需要做第二次消除震荡的步骤。由于$\cos(2\arctan v)=\frac{1-v^2}{1+v^2}$位于$F$表达式的分母中，它和$du=\frac{dv}{1+v^2}$刚好抵消了，变成了$\int_{-\infty}^{\infty} \frac{\sin^2 v}{a+bv^2} dv$的形式，二倍角公式之后可以看出来这就是要求$\frac{1}{a+bx^2}$的傅里叶变换，其封闭形式前文已经介绍过了。如果非要强行赢的话，我觉得AI这样做是在解决数学竞赛问题，偏离了数值计算的初衷了。
第五题：崎岖表面上最低的针
A needle of length 1 rests on the surface defined by the height function $h(x)=0.1x^2+0.1\sin(6x)+0.03\sin(12x)$. What is the lowest possible height of the center of the needle?
“针中点的最低高度”，为什么要优化这么奇怪的东西呢？把它想象成一个物理问题就好了，如果针是均匀的，这个问题就是要求让针的势能最低的位置。首先从图形上理解一下这个曲线是什么样的：
可见这一曲线是相当的崎岖，想直接求出最优解恐怕不太可能。与第三题类似，我决定先去找一个近似解，然后再针对这个近似解的性质研究如何优化它。首先将整个问题离散化，枚举针的全部可能角度，这一问题就变成了
已知有$n$个元素的一维数组$a$，给定斜率$k$，要求对于每一个长度为$m$的区间，找到最小的$b$使得区间内$a_i\le ki+b$。
虽然比第一题还是难一些，但这个就算放在普及组也算是水题了，直接单调队列解决，C++代码。我将它求出的近似解绘制如下：
由此可见这根针下方与曲线相切，而右端顶在曲线上。使用这个近似解计算可得切点的$x$坐标约为$x_0=0.23684$，而右端和曲线交点的$x$坐标约为$r_0=0.90241$。虽然整体来看这一优化问题的目标函数相当崎岖，但在最优解的附近它当然是凸的。因此我们可以直接在$x_0$的附近三分法搜索切点，然后使用$r_0$为初值牛顿迭代计算切线和曲线的交点，从而确定整个线段的位置，并计算中点的高度。C++代码，答案为
\[\boxed{0.0768977459.}\]AI在这道题的表现也相当拉垮。它一开始没有理解”rest on the surface”是什么意思，给出的解位于曲线的下方。我要求针位于曲线上方之后，AI使用了错误的假设：针的两端一定都位于曲线上。我如果只告诉它这个解是错的，它要么就认为这是一个脑筋急转弯，开始揣摩文字的其余意思；要么就是以为是它的解不够精确，需要refine。只有在我点出来它的假设是错误的之后，它才能成功求解。AI最后的做法与我的做法本质上是类似的，它把最优化的目标函数具体写出来了，将其求导，然后把整个区间等分为几千个网格搜索零点，最后在粗筛出的候选网格上用二分法求解高精度答案。从这道题和第三题中可见，AI对于需要图形思维帮助的问题，目前的表现还远不及人类。
第六题：不含42的自然数倒数和
What is $\sum n^{-1}$, where $n$ is restricted to those positive integers whose decimal representation does not contain the substring 42?
与这道题相关的著名结论：十进制表示不含9（或者任意给定子串）的自然数倒数和是收敛的。这是我觉得最难的一道题，我最后想出的解法也是最菜的。对于这种性质很不好、收敛又很慢的数值求和，我死活想不到该怎么办。在精确计算之前，先估计一下大概要算多少项能达到所需的精度吧。在$[10^k,10^{k+1})$之中符合条件的数占比约为$(99/100)^k$，而调和级数的和是$\log(10^{k+1})-\log(10^{k})=\log{10}$。如果我们假设这些数是均匀分布的，求和结果即为$\sum_{k\ge 0}\log 10(99/100)^k$，达到10位精度大概需要$k=3000$。
这个估计结果当然与真实值有很大的偏差。偏差有两个来源：
- 符合条件的数占比并非$(99/100)^k$；
- 符合条件的数不是均匀分布的，例如$\overline{42\dots}$这一个很大的区间对求和都没有贡献。
第一个偏差很好解决，用递推容易求出不含42的数的精确占比。而第二个偏差就不好解决了：这一分布相当不规则。我经过大量尝试之后注意到，如果固定$n$的位数以及其十进制表达的前8位数字，在这个集合中不含42的数字虽然不是均匀分布，但强行视为均匀分布的估算结果和真实值误差小于$10^{-10}$。因此可以枚举位数，再枚举前8位，在被枚举的区间$[L,R]$中使用$\log(R)-\log(L)+\frac{1}{2}(\frac{1}{L}+\frac{1}{R})$作为调和级数的精度足够高的近似（其误差为$O(L^{-2})$，而$L\ge 10^7$，因此近似精度是够的。这个梯形公式本身是平凡的，更高阶的近似可以用Euler–Maclaurin formula做），用递推出来的不含42的数的精确占比乘上调和级数再求和就得到了答案。
然而由于需要枚举的位数高达$3000$，这一做法的用时仍然是不可接受的。当后面填充的位数$j$较高时，区间$[i\times 10^{j},(i+1)\times 10^{j}-1]$上的调和级数部分和完全不需要用梯形公式去做到高精度，而可以直接用$\log(i+1)-\log(i)$近似：这个数是与$j$无关的！实践证明只有$j<3$时才需要梯形公式来求和。因此我们可以直接预处理递推数组的后缀和，前面3项用老办法处理，然后加上后缀和与$\log(i+1)-\log(i)$的乘积即可。C++代码，答案为
\[\boxed{228.4463042.}\]这道题是我被AI虐得最惨的一道，我将此归因为它是所有的问题中公开资料最多的一个，甚至有自己的wiki页面：Kempner series。总之对AI来讲这不过是又一道原题罢了，它直接写出了一个优秀的动态规划，可以把答案计算到任意高的精度，完爆了我的指数级复杂度算法。
在讲解AI是怎么吊打我的之前，我需要介绍一下一种常用的数值计算加速思路：拆贡献。我第一次接触它是在初三时有一道OI题要在区间操作中维护$\log(x)$的求和；这当然是无法维护的，但如果把$\log(x)$泰勒展开成$O(\log \epsilon^{-1})$个$(x-1)^k$的求和，就变得能维护了。当时做出这道题给了我很大的成就感，至今难忘。所谓的“拆贡献”就是利用数值近似，将不可处理的函数$f$分解为可处理的一组基底函数的和。FMM是这样的，Fast Gauss Transform也是如此，而应用范围广泛得多的SVD与低秩近似本质上也是在“拆贡献”。
回到这道题上来，我们需要的求和是$\frac{1}{a\cdot 10^m+b}$，其中$a$和$b$的十进制表示不出现42。该求和本身无法直接处理，可是对于任意的自然数$j$，$\sum_b b^j$是能被递推出来的，因此我们需要用这组基底来表示$\frac{1}{a\cdot 10^m+b}$：泰勒展开即可。具体来讲，由于
\[\frac{1}{a+x}=\sum_{j\ge 0}(-1)^j\frac{x^j}{a^{j+1}},\]$\sum_b b^j$乘上的系数是$\frac{(-1)^j}{10^m}\sum_a \frac{1}{a^{j+1}}$。当$a\cdot 10^m\gg b$时，很少的$j$就足够达到高精度了。让$a$有5位数就可以保证这个“远大于”，而且范围也足够小，可暴力计算$\frac{1}{a^{j+1}}$的求和。我们只需枚举$O(\log \epsilon^{-1})$个$j$，将$a$的求和乘上递推出的$\sum_b b^j$，对全部$j$求和即可。
这个分拆并不是什么新鲜的思想，其实和CF 1578G一模一样。当然，考虑到Trefethen教授能把那种问题放到第一题，我有充分的理由相信这不可能是他的intended solution。但确实该问题的多项式算法在2008年就被人研究过了，见Summing a Curious, Slowly Convergent Series这篇论文，其讨论了不含任意给定子串的自然数的倒数和。文中的做法相当繁琐，被AI给出的优美算法吊打；而且它在重新发明轮子，针对字符串构造了一个“状态之间的有向图”，其实就是KMP的失配函数形成的DFA，但显然论文作者对此一无所知。显然把AI的做法改成在DFA上dp求$\sum_b b^j$就能解决这个一般的情况了。希望AI时代这种隔行如隔山导致的笑话会少一些吧。
第七题：不规则区域的面积
If $f(x,y)=exp(-(y+x^3)^2)$ and $g(x,y)=\frac{1}{32}y^2+e^{\sin y}$, what is the area of the region of the $x$-$y$ plane in which $f>g$?
这道题应该是所有问题中第二简单的了，但隐藏着小小的陷阱。首先我用GeoGebra绘制了这个区域，准备针对区域的形状决定该如何求解。
很好，这个区域的边界看起来相当光滑，而且是星形的，这就非常简单了。我们在星形的中心随便选择一个点为原点建立极坐标系，面积就是$\frac12 \int_{0}^{2\pi}r(\theta)^2d\theta$，数值积分即可；其中的$r(\theta)$表示射线和边界的交点，容易用二分法计算。
真正有趣的问题是这个数值积分该如何计算。我们学习了很多高端的数值积分算法，但对这个问题，其实最高效的算法反而是最朴素的：直接取$n$个等距的$\theta$加起来就是积分的良好近似了。感性理解，这是由于圆周的对称性，任何特殊的采样和权值设置方法都是没有道理的；严谨证明也很简单，考虑三角多项式的正交基$\{e^{ik\theta}\}$，根据单位根反演显然次数小于$n$的三角多项式在$n$次单位根上的和都与积分相同，因此是代数精度最高的$n$节点积分方法。施加变换$z\to \operatorname{Re}(z)$之后这就变成了$[-1,1]$上的Chebyshev多项式与相应的积分公式。这是我认为最优美的数值分析结论。
那么具体的收敛速度是什么呢？对解析函数来说是线性收敛，即$O(\log \epsilon^{-1})$，因此算个10位小数完全不成问题。为什么收敛这么快呢？由于Cauchy积分公式，我们可以把傅里叶级数——即乘上$e^{ik\theta}$的积分“平移”到$\theta+i\sigma$上处理，相当于乘上$e^{ik(\theta+i\sigma)}$再积分，因此傅里叶级数的衰减速度为$e^{-k\sigma}$，而前$n$阶的傅里叶级数部分和的积分是精确的，误差自然是$O(e^{-n\sigma})$，其中$\sigma$的取值取决于奇点离实轴有多远。这个结论在本书中也提到了（虽然是说的区间上的Gauss-Legendre求积，但本质上是同理的），Trefethen教授吐槽道：所有数值分析课程都会教授Gauss-Legendre求积，但从未讲过如此重要的结论，只是因为这需要把解析函数延拓到复平面上来证明，而数值分析课程从不触及复变函数！对此我也深有同感，复分析在数值分析课中被忽视太多了，比如我本科毕设研究的特征值问题，很多新的快速算法都是从围道积分来入手的，由于解析函数在复平面上的优美性质，围道积分总是能带来实数方法无法比拟的收敛速度。
这道题到此为止似乎就完满地解决了，可我在与均匀取点求和的暴力法对拍的时候发现：从图中看，$f>g$的区域完全位于$[-1,2]\times [-3,0]$之中，可是我在该区域外调整暴力搜索的区域边界时，答案却居然随着边界的选择而改变！改变的幅度虽小，但其实远超截断误差所允许的理论上限。仔细一想后我意识到这一现象的原因：我们只绘制了一个$f>g$的区域，可真的只有这一个区域吗？答案是否定的。我让暴力程序输出它找到的其它$f>g$的点，结果发现居然在$(-1.6,4.1)$附近还有另一个$f>g$的区域！
这个区域实在是太小了，GeoGebra在默认的显示区域绘制$f=g$的图时甚至压根没有显示出来，不然我也不会漏掉了。很容易修复，把之前的代码在这个区域再跑一次加起来就行了。C++代码，答案是
\[\boxed{2.115029881.}\]这道题上AI倒是一次性就做对了，略优于我。无法利用图形思考在这道题中对AI反而是个优势，由于这道题中$x$其实只出现了一次，直接将不等式变形就能用$y$表示出$x$的取值范围，转为一维积分就随便算了，不会有漏掉区域的情况。如果把该题中的$f$和$g$设置得更加复杂一些，让任何化简操作都不可行，我相信会和第三题与第五题一样被人类吊打。
第八题：两个相邻正方体的引力
Two adjacent solid unit cubes, each with mass 1, attract each other gravitationally according to Newton’s Law with constant $G=1$. What is the force of attraction between them?
前面的题AI理解有错属于AI的问题，但这题是真的语义不明，我必须补充一下条件：方块是均匀的；求的“引力”是所有质点对的引力之和，不考虑潮汐力之类的高阶效应。
如果方块不是连续的，而是$10^6$个离散点的话，这完全就是标准的FMM。连续和离散并没有本质的区别，两个区域之间的相互作用只需要按照“远场”和“近场”来分类处理。远场之间的作用是低秩的，换句话说每一对点的贡献都能被一个低次多项式良好近似，我们就是要求这个低次多项式的积分，用Gauss-Legendre求积很容易就能做到10位有效数字。而近场的部分递归下去做就好了。
我在实际去做之前还是尝试稍微优化一下。设两个立方体的底面都是$[0,1]^2$，我们要计算立方体$[0,1]^3$对立方体$[0,1]^2\times [0,-1]$对内所有点的引力积分，相当于对所有$(x,y)\in [0,1]^2$求沿着$z$方向的引力积分，再在$[0,1]^2$上积分；由于引力的合力是沿着$z$方向的，我们只需要积分$z$方向的分量，因此一条线段上的积分结果就是两个端点的引力势能之差。于是，我们只需要计算$([0,1]^2,-1)$上的引力势能积分，减去$([0,1]^2,0)$上的引力势能积分，把积分从6维降到了5维。并且，最后1维的积分$\int_0^1 \frac{1}{\sqrt{c+(z-z_0)^2}}dz$是有封闭形式$\operatorname{arsinh}(z_0c^{-1/2})+\operatorname{arsinh}((1-z_0)c^{-1/2})$的，因此只需要四维积分就足够了。其实还能进一步优化，但我当时没想到。
在$([0,1]^2,-1)$上的引力势能积分是远场，可以直接数值积分；而立方体在其底面上的引力势能是近场，需要递归下去。我们把立方体均匀切成上下两层，则上层对底面是远场，可以直接积分；而下层对底面的积分可以看成四个小立方体对底面的积分，加上8个小立方体对底面的相邻方格的积分，加上四个小立方体对底面的对角方格的积分。由于规模为$1/2$的立方体在一半距离处的引力是原问题的$1/4$，而积分面积也是原问题的$1/4$，其答案为原问题的$1/{16}$。如果我们设立方体$[0,1]^3$产生的引力势能在$[0,1]^2$上的积分为$x$，在$[0,1]\times[1,2]$上的积分为$y$，在$[1,2]^2$上的积分为$z$，则根据上面的分析，有
\[x=x/4+y/2+z/4+c_0,\]其中$c_0$是可以用数值积分算出的常数。类似地，有
\[\begin{aligned} y&=y/8+z/8+c_1,\\ z&=z/{16}+c_2. \end{aligned}\]从方程组中可以解出$(x,y,z)$，这题就做完了。C++代码，与第四题一样使用了现成的复合五点Gauss-Legendre求积代码。答案为
\[\boxed{0.9259812606.}\]Trefethen教授是直接对六重积分采取我上面说的分割与递归方法来求解的。其实对于任意的长方体这都是能求出解析解的，而且解法的时间远早于该问题的提出，具体请查看论文The gravitational attraction of a right rectangular prism。
AI对这题的解法和我的各有优劣：它注意到被积函数只和$\mathbf{r_1}-\mathbf{r_2}$有关，因此直接换元$\mathbf{r}=\mathbf{r_1}-\mathbf{r_2}$，降阶为了3重积分。$\mathbf{r}$的分布是什么呢？信号处理里众所周知，box function和自身的卷积是tent function，3维相互独立，因此是$\mathbf{1}-|\mathbf{r}|$。而对于这个3重积分，AI成功意识到了任何一个分量在0附近都会有奇性，因此把三维都分成了多段，越靠近0部分分段越密；在每个小块上用很高阶的Gauss–Legendre来硬做求积达到足够的精度。这种自适应的求积方法对于实际问题非常有效，但对于这种优雅的问题就比利用问题性质的做法要低效不少了。其实我的4重积分也只和$\mathbf{r_1}-\mathbf{r_2}$有关，能用这一手段降到2重积分；这条路再往下走最后能通向极为复杂的解析解，但我对此并无兴趣。
两个微分方程问题
由于我不了解微分方程，没有真正去做，只给出了一些初步设想。
What is the smallest value of $\epsilon>0$ for which the equation $\epsilon u''+u-u^3=0$ with $u(\pm 1)=0$ has exactly five solutions?
为了求解这个问题显然需要一个能把微分方程的数值解计算到10位有效数字的办法，但我对此一无所知，因此就没做。在有了这个求解器之后，这个问题本质上就是对每个$\epsilon$给出了一个$\mathbb{R}\to \mathbb{R}$的连续映射$f_{\epsilon}$，把$u’(-1)$映射到$u(1)$，需要$0$有5个原像。如果真有求解器的话，绘制$\epsilon\mapsto f_\epsilon^{-1}(0)$的图会很有意思，可能类似Logistic map，具有一系列的分歧点。然后从图上就能看出大概什么位置是有5个原像的，二分求解即可。
At what time $t_{\infty}$ does the solution of the equation $u_t=\Delta u+e^u$ on a $3\times 3$ square with zero boundary and initial data blow up to $\infty$?
对ODE我还可以说是一知半解，但是PDE就真的是一无所知了。这个求解器会比前一个问题的复杂非常非常多，甚至感觉都能发一篇论文的程度，需要不仅在平缓变化的区域精度有保障，还能自适应划分网格，在奇点的附近仍然保持极高精度。总之难度全在solver上，有了solver这一问题就是平凡的。
