这是一个关于我如何能够将长时间运行的程序搜索大量整数的运行时间缩短约 95% 的故事。虽然这个特定问题可能听起来很深奥，但这种类型的事情会在许多科学和金融环境中出现，在这些环境中，往往会大量使用多精度算术库。基本的见解无非是放弃 GNU MP 以支持内置数据类型。虽然在必须使用多精度算术库的情况下，GMP 没有很好的替代方案，但当您需要不超过 128 位时，它可能会产生过多的开销。一段时间以来，brian d foy 和我一直在交流他在逐渐扩大的空间中寻找优秀数字的事情。我必须承认，到目前为止，我对寻找 π 的另一个数字或寻找另一个素数之类的事情并不感兴趣。但是，像往常一样，布赖恩让它变得更有趣。如果存在两个正整数 a 和 b 使得： 其中 a 和 b 具有相同的位数 d 且 K 为 10 d/2，则称正整数 n 是一个极好的数。例如，48 是一个极好的数字，因为它等于 8 2 - 4 2。另一方面，24 不是，因为 4 2 - 2 2 是 12。brian 探索了关于极好的数字的各种事实。他还用多种语言编写了程序来搜索这些数字。
这个问题一开始听起来很简单，但随着一个人从 10 位数字毕业到 20 位及以上的数字，它变得越来越难。正如 Mark Jason Dominus 差不多十年前指出的那样，通过不搜索整个 d 位空间，而是通过 d/2 位空间寻找我们可以找到 ab 的 a，可以节省大量工作。这使得连接 ab 和极好的数字。 brian 做了额外的工作来寻找 a 的上限，这样就找不到具有相同位数的候选 b。这些和其他一些技术最大限度地减少了所涉及的工作。但是，人们必须应对所涉及的规模。您现在不想在处理非常大的数字时必须考虑浮点错误，是吗？ Perl with Math::GMP 太慢了，所以 brian 转而使用 C，并使用 GNU MP 编写了一个程序。他在运行 1.8GHz 英特尔酷睿 i5 的 MacBook Air 上获得了以下时序：鉴于这些时序，brian 为他的程序配备了一些细节，以便该程序可以提供定期进度报告，通过报告其取得的进展干净利落地终止，以及从给定点重新开始搜索等。
在这一点上，虽然通过在具有数十个内核的机器上运行许多实例似乎可以发现所有 30 位的优秀数字，但很明显，即使在 32 位空间中制造一个小凹痕也需要花费大量时间核心，或重大发展。经过我们的交谈，我的兴趣又被重新点燃了。起初，我尽我最大的努力回忆起我大约 20 年前研究的抽象代数和代数拓扑的一些内容，看看我是否能找到一种方法，甚至略微改进搜索空间的限制。当我无处可去时，一个想法在我的脑海中浮现：存储一个 16 位整数需要多少位？ Easy peasy：16×log 2 10 大约是 53.15，所以超过 53 位就足以存储 32 位整数的每一半。存储一个 18 位整数只需要不到 60 位，因此，如果我可以只进行整数数学运算而不是依赖 GMP，那么最多 36 位数字可能会有一些显着的速度提升。在那个 sqrt 中有两个 64 位数字相乘的两个实例。我们需要 128 位来保存中间值，我们最终将一个 128 位整数转换为双精度值。这避免了处理 128 位中间值。根据定义，a 的范围从 10 (d/2 - 1) 到 10 d/2。因此，平方根内除法的分子和分母就像大小一样，这应该可以减少将 64 位整数转换为两倍并除以它们的浮点误差，但我必须承认，我还没有坐下来计算误差范围.一旦我们有了给定 a 的候选 b，我们就会返回并检查是否满足优秀数字的定义，就像 Brian 已经在他的基于 GMP 的程序中所做的那样。当然，这实际上涉及将两个 64 位整数相乘几次。
当我想到这一点时，我正在使用 Windows 旧车，所以我的第一次尝试是使用 Visual Studio 2015 附带的 C 编译器提供的 128 个整数运算编写的。请注意，该编译器不支持 128 位整数类型，但是它确实提供了一个 UnsignedMultiply128 函数，不管它有多麻烦： #include <math.h> #include <stdint.h> #include <stdio.h> #include <stdlib.h> #include <windows.h> const uint64_t powers_of_10 [] = {1，10，100，1000，10000，100000，百万，千万，亿，十亿，百亿I64，千亿I64，万亿I64，十万亿I64，百万亿I64，千万亿I64，10000000000000000 I64，100000000000000000 I64，1000000000000000000 i64, }; int main( int argc, char *argv[ 1]) { int d, k; uint64_t K，开始，前面，后面，last_digit，计数= 0； uint64_t lhs[2]、rhs[2]、frontsq[2]； if (argc < 2) { fputs( "需要位数", stderr);退出（1）； } d = atoi(argv[1]);如果 (!d) { d = 2; } if (d % 2) { d *= 2; } k = d/ 2; if (k >= sizeof(powers_of_10)/ sizeof(powers_of_10[ 0])) { fputs("数字太多", stderr);退出（1）； K = powers_of_10[ k ];开始 = powers_of_10[k - 1]; for (front = start; front < 7 * start; front += 1) { last_digit = (front % 10); if ( ((last_digit != 0) && (last_digit != 4) && (last_digit != 6)) { continue; } back = ( uint64_t) ( 1.0 + front * sqrt( 1 + ((double) K) / front));如果（返回 >= K）{ 中断； } lhs[ 0] = lhs[ 1] = rhs[ 0] = rhs[ 1] = frontsq[ 0] = frontsq[ 1] = 0; lhs[0] = UnsignedMultiply128(back, back - 1, lhs + 1); rhs[0] = UnsignedMultiply128(front, K, rhs + 1); frontsq[0] = UnsignedMultiply128(front,front,frontsq + 1); rhs[0]+=frontsq[0]； rhs[1]+=frontsq[1]； if (rhs[ 0] < frontsq[ 0]) { rhs[ 1] += 1; } if ((lhs[ 1] == rhs[ 1]) && (lhs[ 0] == rhs[ 0])) { count += 1; printf( "%I64d%I64d \n ", 前面, 后面); } } printf( "%I64d 具有 %d 个数字的优秀数字 \n ", count, d);返回0；请注意，与 brian 计算的相当严格的界限相比，该程序采用了非常粗略的空间缩减策略，因此它会浪费时间在杂草中徘徊。它还检查连续的 a 值，而不是使用跳转表等。我首先用八位数字检查它。它的输出与 brian 找到的数字匹配，所以我跳到了 14 位数字：这个东西上的 CPU 是一个相当陈旧的 Intel Core 2 Duo T7600 @ 2.33 GHz 与 1.8 GHz i5 相比，brian 用他的 GMP 估计时间程序。然而，他花了大约一秒钟的时间穿过 14 位数字的空间，而这个笨拙的小东西只花了不到 1/5 的时间。
然后，与基于 GMP 的程序花费的 3.5 小时相比，我在大约 17 分 5 秒内完成了 22 位数字空间。最后试了一下24位空格，只用了2小时49分钟。相比之下，基于 GMP 的程序需要 36 小时。我欣喜若狂，就这个好消息联系了布赖恩。这意味着在我们对这些数字的代数结构的理解上没有巨大突破的情况下，超越 30 位空间最终是可行的。当然，考虑到长时间运行的任务是在运行 Linux 的多核计算机上执行的，没有太多理由坚持使用 Visual Studio 特定功能。幸运的是，gcc 实际上提供了一个内置的 __int128 类型，它使一切变得不那么混乱。我改编了 brian 的基于 gcc 的程序以仅使用基本类型，然后将其关闭，期待更出色的结果。这是乘法例程：excellent_full_tmultiply_halves(constexcellent_half_t x,constexcellent_half_t y) {excellent_full_t z = ((excellent_full_t) x) * ((excellent_full_t) y); return z;} voidcheck_excellent(excellent_half_t a,excellent_half_t K) {excellent_half_t b = 1.0 + a * sqrt(1 + ((excellent_float_t) K)/ a); excellent_full_t lhs = multiply_halves(b, b - 1); excellent_full_t rhs = multiply_halves(a, K) + multiply_halves(a, a); if ( lhs == rhs ) { print_excellent_number(a, b); } return;} 然后我意识到发生了什么：brian 的程序安装了信号处理程序以提供按需统计数据和响应 CTRL-C 的优雅退出。信号处理程序设置 volatile sig_atomic_t 类型的全局标志，以向内循环指示已接收到信号。基于 GMP 的程序足够慢，以至于在每次迭代时检查这些标志不会产生明显的减速。
然而，使用内置函数，在相同的硬件上，我的程序能够处理更多数量级的候选数字，因此每秒检查的数量也激增。仅供参考，检查 volatile sig_atomic_t 变量的值是一项昂贵的操作。因此，我调整了程序，使其仅在可配置的时间间隔内检查信号。该程序使用它能够执行的每秒迭代次数来调整它检查标志的频率。这些都不是很准确，但也没有那么重要。我们只希望程序能够以合理的方式响应信号，而不会过多地分散其主要任务的注意力。进一步的优化包括使用 -O2 -march=native -ffinite-math-only -fno-math-errno 构建以提供更多性能。 FWIW，其中一些选项在 OSX 上使用 clang 似乎不可用，并且当使用 clang 构建时，GMP 版本似乎没有通过 -O2 和 -march=native 改进。事实上，在 ArchLinux VM 中使用 gcc 构建的相同程序在同一 Mac 上的 VM 中运行速度比使用 clang 构建的本机版本更快。我想在同等条件下比较基于 GMP 的程序与使用内置程序的程序。为此，我使用了配备 Core i3 Broadwell CPU @2.1Ghz 的华硕笔记本电脑。我还没有机会在这台机器上安装 ArchLinux，所以必须使用 Cygwin 的 gcc 4.9.3。
公平地说，为了进行比较，我从两个程序的内部循环中删除了信号检查代码。它对我的程序没有太大影响，因为信号检查默认为每两秒一次，但它使基于 GMP 的程序运行得更快一点。这是并排比较：