全靠数学技巧，在60毫秒内计算十亿级阶乘？ 是的，你没看错。 一位名为Adamant的网友就展示了他如何将十亿级阶乘的计算时间，从数秒压缩到惊人的61毫秒！整个过程无需预计算，也无需FFT。 这个挑战是这样的：给定 𝑛1,…,𝑛𝑡 ，其中 𝑡≤105 。对于每个 𝑖 ，求 𝑛𝑖!mod 𝑀 ，其中 𝑀=998244353。 如果我们暴力硬解，比如【图1】，计算需要耗时3745毫秒。 但如果我们用上下面这些技巧，就能实现60倍提升，将时间压缩至61毫秒。不卖关子，直接看： 第一步：使用威尔逊定理 威尔逊定理有个变式：𝑛!(𝑝−1−𝑛)!≡(−1)ⁿ(mod 𝑝)。 如果选择𝑛和𝑝−1−𝑛中较小的那个，就只需要计算到𝑝/2，而不需要遍历所有𝑝个数。 这样做可以将运行时间减少到1889 毫秒 ，预计可实现2倍的提升。 第二步：跳过偶数阶乘 我们知道，2⋅4⋅…⋅2n=(2n)!!=2nn!。这意味着我们可以不断提取偶数因子，直到只剩下奇数的乘积，再乘以一个很大的2的幂次。 设𝑓(𝑛)为所有不超过𝑛的奇数的乘积。阶乘可以重写为𝑛!=∏ₖ₌₁^∞ 2^⌊𝑛/2ᵏ⌋ 𝑓(⌊𝑛/2ᵏ⌋)。 实施这一改变后，运行时间降至998毫秒，又获得了近2倍的提速。 第三步：利用指令级并行 在动用像SIMD和向量化这样的大杀器之前，我们还有招——指令级并行。 现代CPU能同时处理多条指令。虽然阶乘计算中，每一步都依赖上一步的结果，但我们可以并行处理多个独立的计算块。 例如，通过将计算分成8个独立的块并行处理，运行时间可以进一步降至287毫秒，提速3.5倍。 第四步：向量化的蒙哥马利积 这一步是将计算块细分成更小的子块，并利用向量化技术在每个子块中累加结果。【图2】 由于蒙哥马利乘法将(𝑎,𝑏)转换为𝑎𝑏⋅2⁻³²，我们可以跟踪乘以的2的总幂次，在最后消除这些累积的2⁻³²，这样就不会在主循环中干扰。 这样做之后，时间降至119毫秒，又获得了2.4倍的提速。 第五步：快速逆元、幂次和输入 在之前的优化中，对特别大的𝑥计算2ˣ mod 𝑝和为某些结果求逆元会耗费大量时间。通常它们需要O(logp)的时间，但我们可以将其优化到O(1)。 具体做法是： 预先计算2的幂次表，用于快速求2x(mod p)。【图3】 对于多个数的逆元，我们可以在O(n+logp)时间内批量计算，远快于O(nlogp)。【图4】 结合这些优化和更快的输入/输出（IO），运行时间进一步降至99毫秒，提速1.2倍。 第六步：回退到常规阶乘计算 将每个n!的计算都转化为“奇数阶乘”的查询会带来一些存储和乘法累加的压力。 更好的方法是只在最大的数上使用奇数阶乘的优化，当n足够小时，就回退到更常规的阶乘计算方法。 加上这个优化后，速度就达到了64毫秒，实现了最后的1.5倍提速。 最后，在使用一些小优化，就可以把时间推到61毫秒，在此就不多赘述了。 博客原文：