如果已知非常大的整数 N 的因式分解,则可以快速计算浮点 1/N

Art*_*oul 12 c++ python algorithm math performance

如果我有知道分解的整数,那么计算浮点数N最快(最有效)的方法是什么?1/N为此应使用大浮点(或整数)算术。

我想在 C++ 中执行此操作(或在 Python 中进行实验运行)。

我的N非常巨大,大小为千兆/太比特。生成的浮点N还应该具有很高的精度,大约与初始的位大小相同N

需要精确的浮点值,这意味着如果我请求浮点精度Log2(N)位,那么至少95%结果的所有前导位应该是精确的(所有位与理想值中的位相同)。

当然,4^Ceil(Log2(N)) / N如果它有助于和/或简化任务,则可以计算整数除法,而不是浮点计算。对我来说,这两个任务(整数和浮点数)本质上是相同的,因为整数表示可以转换为浮点数,反之亦然。

一个重要的注意事项是 的因式分解N只有很小的质因数,它们的大小都是 32 位(当然最多可能是 64 位)。

我想知道如果因式分解N以及因数很小这一事实是否可以帮助以某种方式解决任务?

当然,我没有首先实现我自己的部门,而是尝试使用高度优化的GMP库来完成此任务,但它(据我所知)并没有使用N已经因式分解的事实。

任何人都可以建议我是否要为此实现自己的函数,只是为了通过实验弄清楚它是否比 GMP 更快,那么我应该使用哪种算法?

我发现这里可以使用 3 种算法 1)长除法,这是一种学校级算法。2)巴雷特还原。3)蒙哥马利还原

其实我不知道其他算法。你能建议其他吗?仅当相同素因子重复多次时,Barrett 和 Montgomery 约简才有帮助,否则单次除法不值得 Barrett 和 Montgomery 所需的预计算。

此外,Barrett/Montgomery 归约仍然需要一次性计算4^Ceil(Log2(N)) / PrimeDivisor。所以它们不会让你免于进行长除算法。

对于长除法算法,我肯定会使用它2^64作为基础而不是基础10(如在学校)。

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

当然,在长除法中我可以使用 Montgomery 和 Barrett,因为在每一步都需要短除法(128 位整数除以 64 位整数)(如果它能提供任何提升)。

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

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

Jér*_*ard 5

大概的概念

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

问题是,1/a许多因素的计算a都可以比计算更快1/N。的确:

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

然而,有一个问题:该方法需要大量乘法,因为存在很多因素,尽管它们大部分可以并行计算,但这仍然是一个相当昂贵的计算(尤其是最后一个乘法,很难在平行线)。因此,我不确定该方法比 GMP 中实施的非常优化的方法快得多,但它确实值得一试。


笔记

最好将相同大小的数字相乘,甚至将第一个数字与最小的十进制/二进制扩展周期相乘(为了最小化计算倒数的符号表示的大小,从而导致更小的内存占用和更快的速度)乘法)。

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

您可以使用平方求幂算法来有效计算很大p^e时的倒数。当很大时e,这样的项经常出现在素数分解的指数形式中。NN

如果十进制/二进制扩展的周期大于最终的浮点数,您可以截断它。但这可能会影响该方法的准确性。我认为至少需要一些额外的数字才能获得准确的结果(特别是由于应用了许多乘法)。

对于涉及大量符号表示的最后乘法,您可以生成大浮点数并将其提供给 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
Run Code Online (Sandbox Code Playgroud)

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

  • “1/a”的周期长度最大为“a-1”,例如 375292573(质数)的倒数的二进制周期长度为 375292572。所以这可能行不通。 (2认同)
  • @Arty 1/a 的周期长度是 [2 mod a 的乘法阶数](https://en.wikipedia.org/wiki/Multiplicative_order)(2 是基数;对于基数 10 的分数,您需要阶数10 模 a)。[此](https://rosettacode.org/wiki/Multiplicative_order) 是多种语言的乘法阶计算器(不保证很快)。 (2认同)