如果已知一个非常大的整数N的因数分解,如何快速计算浮点数1/N?

12
如果我已知正整数 N 的因数分解,那么计算浮点数 1/N 最快(最高效)的方法是什么?需要使用大浮点数(或整数)算术。
我想在 C++ 中完成这个任务(或者在 Python 中进行实验运行)。
我的 N 非常庞大,大小为千兆/兆位。得到的浮点数 N 也应具有巨大的精度,约为初始 N 的位数。
需要精确的浮点值,这意味着如果我请求浮点精度为 Log2(N) 位,则结果的至少 95% 的前导位应该是精确的(与理想值相同)。
当然,如果有帮助和/或简化任务,可以计算整数除法 4^Ceil(Log2(N)) / N,而不是浮点计算。对于我来说,这两个任务(整数和浮点)本质上是相同的,因为整数表示可转换为浮点表示,反之亦然。
一个重要的注意事项是,N 的因数分解只有小的素因子,它们都是 32 位大小(最多可能是 64 位)。
我想知道是否拥有 N 的因数分解以及因子很小的事实,是否可以在解决任务时有所帮助?当然,我会首先尝试使用高度优化的GMP库来实现这个任务,而不是自己实现除法,但是(据我所知),它没有利用已经因式分解的事实。
如果我要实现自己的函数来解决这个问题,并想通过实验弄清楚它是否比GMP更快,那么有人可以建议我使用什么样的算法吗?
我发现这里有3种算法可以使用:1)长除法,这是一种学校级别的算法。2)Barrett约减。3)Montgomery约减
实际上,我不知道还有其他算法。你能推荐另外的算法吗?无论是Barrett还是Montgomery约减都只有在相同的质因子重复多次时才有帮助,否则单次除法不值得为Barrett和Montgomery预计算所需的工作量。
此外,Barrett/Montgomery约减仍然需要一次性计算4^Ceil(Log2(N)) / PrimeDivisor。所以它们并不能避免使用长除法算法。
肯定对于长除法算法,我将使用基数2^64而不是学校里的基数10

我已经使用长除法和其他整数算法以及椭圆曲线算法实现了自己的实验性库。目前它是通用的,因此不比GMP更快。现在我需要一种特殊的除法算法,它至少可以比GMP快几倍。

当然,在长除法中,如果有任何提升,我可以使用蒙哥马利和巴雷特,因为每一步都需要短除法(128位整数除以64位整数)。

此外,在长除法的每个步骤中,我可以使用快速傅里叶变换数论变换进行乘法运算。

以上是我知道的唯一优化方法。是否还有其他可能的优化方式?也许FFT可以用于直接进行除法(而不仅仅是用于乘法)?


我建议您首先让算法产生正确的结果,无论运行速度有多慢。然后在其上运行分析器,看看需要关注哪些缓慢的部分。没有看到实际代码很难给出性能方面的建议。 - Eyal K.
3
@EyalK。目前的问题不是如何优化我已经有的代码。我已经有了通用代码。但是它并不比GMP的通用变体更快,这是显而易见的。优化当前算法的速度没有意义,因为通用解决方案不会比GMP更快。这里的问题是如何利用分解可用性来实现真正快速的算法。同时,也许有人可以告诉我是否有基于FFT的除法技术可用。我不需要只比GMP快2倍的解决方案,如果分解有帮助,我需要一些特殊的算法,可以比GMP快10-20倍。 - Arty
如果你需要无数位小数,我怀疑任何算法都不会比简单的长除法更好(高效、快速)。对我来说,拥有因式分解并不是一种改进,因为它不能避免计算出如此巨大数量的数字和合并部分结果。 - Ripi2
2
例如,对于乘法,存在FFT(傅里叶)算法,可以实现O(N Log(N))的时间复杂度,这比学校级别的O(N^2)算法要好得多,如果你有太比特位数的数字,那么这是一个非常巨大的改进。但是FFT只对大数字快速。也许类似地,对于长除法,存在真正快速的算法,比朴素的长除法更快,但仅适用于非常大的数字。我不知道是否存在这样的算法,但不久之前我也不知道FFT的存在。所以等待数学/算法大师的回复! - Arty
1
你可以在https://math.stackexchange.com/找到更好的数学专家。 - Ripi2
显示剩余6条评论
1个回答

5

总体思路

可以将 N 分解成小的整数会有所帮助。事实上,这意味着 N = a1 * a2 * ... * an。因此,1/N = 1/(a1 * a2 * ... * an) = (1/a1) * (1/a2) * ... * (1/an)

问题是对于许多因子 a 计算 1/a 可能比计算 1/N 更快。实际上:

  • 十进制/二进制展开周期应该远小于1/N,因为log2(a)log2(N)小得多,这意味着1/a的精确表示可以以非常紧凑和快速的方式存储;
  • 因子的倒数可以独立计算,因此可以并行计算;
  • 当一个因子在N的分解中被使用多次时,它昂贵的倒数只需计算一次;
  • 乘法运算应该比计算倒数或除法更快;
  • 最终基于乘法的规约可以使用成对规约策略并行完成。

然而,有一个问题:该方法需要大量的乘法,因为有很多因子,虽然它们可以大部分并行计算,但这仍然是一项相当昂贵的计算(特别是最后一个乘法,难以并行执行)。因此,我不确定该方法是否比GMP中实现的高度优化的方法快得多,但肯定值得一试。


注意事项

最好将相同大小的数字相乘,甚至可以将具有最小十进制/二进制扩展周期的第一个数字相乘(以便最小化计算倒数的符号表示的大小,从而减少内存占用和加快乘法速度)。

32位中大约有2亿个质数,但大多数N的因子应该适合16位(因为较小的因子更频繁),并且只有大约6500个质数适合16位,因此如果您计划计算不同N的多个倒数,则可以预先计算它们的倒数。

e很大时,您可以使用平方取幂算法有效地计算p^e的倒数。当N很大时,在其质因数分解的指数形式中经常出现这样的术语。

你可以截断小数/二进制展开的周期,如果它比最终浮点数大。但这可能会影响方法的准确性。我认为至少需要一些额外的数字才能得到精确的结果(特别是由于应用了许多乘法)。对于涉及巨大符号表示的最后一次乘法,你可以生成大的浮点数并将它们提供给GMP,以便它可以使用非常优化的实现(通常基于快速傅立叶变换)。

示例

这里有一个示例,使用小的N = 307230,为了清晰起见采用十进制表示法(小数部分的循环节已经被标出):

1/307230 = 0.0000032548904729355857175406047586498714318263190443641571461120333300784428603977476157927285746834619015070142889691761872213
              ------------------------------------------------------------------------------------------------------------------------------

307230 = 2 * 3 * 5 * (7^2) * 11 * 19

1/2 = 0.50
         -
1/3 = 0.3
        -
1/5 = 0.20
         -
1/7 = 0.14285714
          ------
1/11 = 0.090
          --
1/19 = 0.0526315789473684210
          ------------------

k1 = ((1/2) * (1/3)) * ((1/5) * (1/11)) = 0.0030
                                          --

k2 = k1 * (1/19) = 0.000159489633173843700
                        ------------------

k3 = (1/7)^2 = 0.0204081632653061224489795918367346938775510
                  ------------------------------------------

1/N = k2 * k3

对于最终的乘法,您可以使用精确算法或使用GMP进行快速近似计算。


1
我是其中一个点赞者,但我怀疑这将输给mpn_invertappr,因为中间结果需要精度要求。 - David Eisenstat
2
1/a的周期长度高达a-1,例如375292573的倒数(一个素数)的二进制周期长度为375292572。因此这可能行不通。 - n. m.
@n.1.8e9-我的份在哪里呢m。你知道有没有一种保证快速找到周期长度的确切方法吗?因为如果我逐位计算数字,得到类似于0.123456456456这样的结果,我不能保证456已经是一个周期,因为后面可能还会有其他数字。除了像你说的那样检查所有的a-1位以获取1/a之外,我还能用什么方法找出周期呢? - Arty
2
@Arty 1/a的周期长度是2模a的乘法阶(其中2为底数;对于10进制分数,需要计算10模a的阶)。这个是一个多语言乘法阶计算器(不能保证速度快)。 - n. m.
@n.1.8e9-我的股份在哪里吗?你能回答另一个问题吗 - 如果我预先计算出一个周期内的1/prime,并乘以这些倒数,是否保证我得到精确的最终答案?也就是说,如果我将其与使用多个周期计算的相同倒数的乘积进行比较,是否会得到相同的最终结果?我只是怀疑单个周期是否足以获得理想的无误结果。 - Arty
显示剩余4条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接