一种计算数学常数 e 的有效方法

Posted 2023-02-22

【中文标题】一种计算数学常数 e 的有效方法

【英文标题】：An efficient way to compute mathematical constant e

【发布时间】：2010-06-12 10:15:52

【问题描述】：

由于许多除法运算，常数 e 作为无限级数之和的标准表示对于计算非常低效。那么有没有其他方法可以有效地计算常数？

谢谢！

编辑

在关注了你的一些链接之后，我相信效率来自一种我不熟悉的称为二进制拆分的技术（虽然表示仍然提到系列）。如果有人熟悉它，请随时贡献。

【问题讨论】：

你为什么要计算它而不是使用常量表达式？

听起来您要不断地计算它，而不是只计算一次以达到所需的精度。

因子增长非常快。泰勒并非效率低下。

然后尝试计算 exp(2) - Taylor 仅对特殊情况 exp(1) 有效，因为 pow(1,N) == 1。

@MSalters 但他不需要 exp(2)

【参考方案1】：

由于不可能计算出“e”的每个位，因此您必须选择一个停止点。

双精度：16 位小数

对于实际应用，“尽可能接近 'e' 真实值的 64 位双精度浮点值——大约 16 个十进制数字”已经绰绰有余了。

正如 KennyTM 所说，该值已经在数学库中为您预先计算好了。 如果你想自己计算，正如 Hans Passant 指出的那样，阶乘已经增长得非常快。 该系列中的前 22 项对于计算到该精度来说已经过大了——如果将结果存储在 64 位双精度浮点变量中，则从该系列中添加更多项不会改变结果。 我认为你眨眼所花的时间比你的电脑做 22 次除法的时间要长。所以我认为没有任何理由进一步优化它。

千、百万或数十亿十进制数字

正如 Matthieu M. 所指出的，这个值已经计算过了，你可以从 Yee 的网站下载。

如果您想自己计算，那么标准的双精度浮点数将无法容纳这么多数字。 你需要一个“bignum”库。 与往常一样，您既可以使用现有的众多免费 bignum 库之一，也可以通过构建自己的另一个具有自己特殊特性的 bignum 库来重新发明***。

结果——一长串数字——并不是非常有用，但计算它的程序有时被用作测试“bignum”库软件的性能和准确性的基准，以及检查稳定性的压力测试和新机器硬件的冷却能力。

一页非常简短地描述了the algorithms Yee uses to calculate mathematical constants。

***"binary splitting" article 更详细。 我认为您正在寻找的部分是数字表示： 而不是在内部将所有数字存储为小数点（或二进制点）前后的一长串数字， Yee 将每个项和每个部分和存储为一个有理数——作为两个整数，每个整数都是一长串数字。 例如，假设其中一个工作 CPU 被分配了部分和，

... 1/4! + 1/5! + 1/6! + ... .

而不是先对每个项进行除法，然后相加，然后将一个百万位数的定点结果返回给管理器 CPU：

// extended to a million digits
1/24 + 1/120 + 1/720 => 0.0416666 + 0.0083333 + 0.00138888

那个CPU可以先用有理算术把系列里的所有项加在一起，然后把有理结果返回给管理器CPU：两个整数，每个可能几百位数：

// faster
1/24 + 1/120 + 1/720 => 1/24 + 840/86400 => 106560/2073600

以这种方式将数千个项相加后，管理器 CPU 在最后进行唯一除法，得到小数点后的小数位数。

记得避开PrematureOptimization，以及 总是ProfileBeforeOptimizing。

【讨论】：

这正是我正在寻找的（关于二进制拆分的部分）。干杯！

@Cinnamon 我在搜索从 SO 到我的网站的链接时遇到了这个问题。在问了这个问题之后，我加入了 SO 很长一段时间，我想我会发表评论，确认二进制拆分方法实际上是关键算法。它降低了对从多项式到准线性的级数求和的计算复杂度。

【参考方案2】：

如果您使用的是double 或float，则math.h 中已经有一个M_E 常量。

#define M_E         2.71828182845904523536028747135266250   /* e */

http://en.wikipedia.org/wiki/Representations_of_e#As_an_infinite_series中还有e的其他表示；都会涉及到分裂。

【讨论】：

我想问题的重点是找到一种有效的算法来计算 e（如果你愿意，只是为了这样做），所以使用预定义的常量不是解决方案。否则请参阅我的答案中的链接以获得您需要的所有精度......

【参考方案3】：

我不知道任何比级数的泰勒展开式“更快”的计算，即：

e = 1/0！ + 1/1！ + 1/2！ + ...

或

1/e = 1/0！ - 1/1！ + 1/2！ - 1/3！ + ...

考虑到这些是由计算 the first 500 billion digits of e 的 A. Yee 使用的，我想没有太多优化可做（或者更好，可以优化，但还没有人找到方法，AFAIK）

编辑

一个非常粗略的实现

#include <iostream>
#include <iomanip>

using namespace std;

double gete(int nsteps)

  // Let's skip the first two terms
  double res = 2.0;
  double fact = 1;

  for (int i=2; i<nsteps; i++)
  
    fact *= i;
    res += 1/fact;
  

  return res;


int main()

  cout << setprecision(50) << gete(10) << endl;
  cout << setprecision(50) << gete(50) << endl;

输出

2.71828152557319224769116772222332656383514404296875
2.71828182845904553488480814849026501178741455078125

【讨论】：

请注意，为了计算 x!，您必须计算 (x-1)！第一的。所以这种表示非常适合动态规划。

@Dave：是的，我想您可以存储前一个分母的结果并将其重用于下一个总和。查看更新的答案

setprecision(50)? double 只能精确到 16 位。

@KennyTM: 是的，只是为了放一个大数字，我不记得double 的精确度了：P

@KennyTM：步数与每一步计算的有效小数位数根本不对应。当然 50 在这里可能是多余的，但估计每次迭代产生的有效小数位数似乎比计算迭代本身要困难得多。

【参考方案4】：

This page 对不同的计算方法进行了很好的总结。

这是一个来自 Xavier Gourdon 的小型 C 程序，用于在您的计算机上计算 e 的 9000 个十进制数字。对于 π 和其他一些通过超几何级数定义的常数，也存在类似的程序。

[来自https://codereview.stackexchange.com/a/33019的脱高尔夫球版本]

#include <stdio.h>
int main() 
      int N = 9009, a[9009], x;
      for (int n = N - 1; n > 0; --n) 
          a[n] = 1;
      
      a[1] = 2;
      while (N > 9) 
          int n = N--;
          while (--n) 
              a[n] = x % n;
              x = 10 * a[n-1] + x/n;
          
          printf("%d", x);
      
      return 0;
  

这个程序[当代码打高尔夫球时]有 117 个字符。它可以更改为计算更多位数（将值 9009 更改为更多）并更快（将常数 10 更改为 10 的另一个幂和 printf 命令）。一个不太明显的问题是找到使用的算法。

【讨论】：

我在codereview.stackexchange.com/a/33019/31087 的引文中进行了编辑，因为（现在）强烈不鼓励仅链接的答案

【参考方案5】：

我在 question regarding computing e by its definition via Taylor series 的 CodeReviews 上给出了这个答案（因此，其他方法不是一个选项）。 cmets中建议使用此处的交叉帖子。我已经删除了与其他主题相关的评论；有兴趣进一步解释的可以查看原帖。

---

C中的解决方案（应该很容易适应C++）：

#include <stdio.h>
#include <math.h>

int main ()

    long double n = 0, f = 1;
    int i;
    for (i = 28; i >= 1; i--) 
        f *= i;  // f = 28*27*...*i = 28! / (i-1)!
        n += f;  // n = 28 + 28*27 + ... + 28! / (i-1)!
      // n = 28! * (1/0! + 1/1! + ... + 1/28!), f = 28!
    n /= f;
    printf("%.64llf\n", n);
    printf("%.64llf\n", expl(1));
    printf("%llg\n", n - expl(1));
    printf("%d\n", n == expl(1));

输出：

2.7182818284590452354281681079939403389289509505033493041992187500
2.7182818284590452354281681079939403389289509505033493041992187500
0
1

有两点很重要：

此代码不计算 1, 1*2, 1*2*3,... 这是 O(n^2)，但一次计算 1*2*3*... （这是 O(n)）。

它从较小的数字开始。如果我们尝试计算

1/1 + 1/2 + 1/6 + ... + 1/20！

并尝试将其添加到 1/21！，我们将添加

1/21！ = 1/51090942171709440000 = 2E-20，

到 2.something，这对结果没有影响（double 包含大约 16 个有效数字）。这个效果叫做underflow。

但是，当我们从这些数字开始时，即如果我们计算 1/32!+1/31!+... 它们都会产生一些影响。

此解决方案似乎与 C 在我的 64 位机器上使用其 expl 函数计算的内容一致，使用 gcc 4.7.2 20120921 编译。

【讨论】：

@VerdanSego 不应该用 n=1 来初始化 n 的总和。这是最小公分母产生的结果。例如3 次迭代 1/0!+1/1!+1/2!+1/3!=(3!+3!+3+1)/3!。只是为了清楚起见，这并不是说它有很大的数值差异。

@AlexanderCska，我认为你是对的。循环有 28 个步骤，但应该有 29 个加数，加数为 0！ = 1 缺失。我需要更彻底地检查，但我现在没有时间。随意编辑。

【参考方案6】：

您可能会获得一些效率。由于每一项都涉及下一个阶乘，因此通过记住阶乘的最后一个值可以获得一些效率。

e = 1 + 1/1! + 1/2! + 1/3! ...  

展开方程：

e = 1 + 1/(1 * 1) + 1/(1 * 1 * 2) + 1/(1 * 2 * 3) ...

不是计算每个阶乘，而是将分母乘以下一个增量。因此，将分母保持为变量并将其相乘会产生一些优化。

【讨论】：

确实，this answer 基本上实现了这一点（加上在1/28! 之后删除任何内容并向后迭代以获得更好的性能）

【参考方案7】：

如果您可以接受最多七位数的近似值，请使用

3-sqrt(5/63)
2.7182819

如果你想要确切的值：

e = (-1)^(1/(j*pi))

其中 j 是虚数单位，pi 是众所周知的数学常数（欧拉恒等式）

【参考方案8】：

有几种“spigot”算法可以无限制地按顺序计算数字。这很有用，因为您可以通过恒定数量的基本算术运算简单地计算“下一个”数字，而无需事先定义您希望产生多少位。

这些应用一系列连续的转换，以便下一个数字出现在 1 的位置，因此它们不受浮点舍入误差的影响。效率很高，因为这些变换可以表示为矩阵乘法，而矩阵乘法可以简化为整数加法和乘法。

简而言之，泰勒级数展开式

e = 1/0! + 1/1! + 1/2! + 1/3! ... + 1/n!

可以通过分解阶乘的小数部分来重写（请注意，为了使系列规则我们将1移到左侧）：

(e - 1) = 1 + (1/2)*(1 + (1/3)*(1 + (1/4)...))

我们可以这样定义一系列函数 f1(x) ... fn(x)：

f1(x) = 1 + (1/2)x
f2(x) = 1 + (1/3)x
f3(x) = 1 + (1/4)x
...

e 的值是从所有这些函数的组合中找到的：

(e-1) = f1(f2(f3(...fn(x))))

我们可以观察到每个函数中 x 的值由下一个函数确定，并且这些值中的每一个都在 [1,2] 范围内 - 也就是说，对于这些函数中的任何一个， x 将是1 <= x <= 2

既然是这种情况，我们可以分别使用 x 的值 1 和 2 来设置 e 的下限和上限：

lower(e-1) = f1(1) = 1 + (1/2)*1 = 3/2 = 1.5
upper(e-1) = f1(2) = 1 + (1/2)*2 = 2

我们可以通过组合上面定义的函数来提高精度，当一个数字在上下界匹配时，我们知道我们计算的 e 值精确到那个数字：

lower(e-1) = f1(f2(f3(1))) = 1 + (1/2)*(1 + (1/3)*(1 + (1/4)*1)) = 41/24 = 1.708333
upper(e-1) = f1(f2(f3(2))) = 1 + (1/2)*(1 + (1/3)*(1 + (1/4)*2)) = 7/4 = 1.75

由于第 1 位和第 10 位数字匹配，我们可以说 (e-1) 的 10 位精度的近似值是 1.7。当第一个数字在上限和下限之间匹配时，我们将其减去然后乘以 10 - 这样，所讨论的数字始终位于浮点精度较高的 1 的位置。

真正的优化来自线性代数中将线性函数描述为变换矩阵的技术。组合函数映射到矩阵乘法，因此所有这些嵌套函数都可以简化为简单的整数乘法和加法。数字相减乘以10的过程也构成线性变换，因此也可以通过矩阵乘法来完成。

方法的另一种解释： http://www.hulver.com/scoop/story/2004/7/22/153549/352

描述算法的论文： http://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf

通过矩阵算术执行线性变换的快速介绍： https://people.math.gatech.edu/~cain/notes/cal6.pdf

注意，该算法利用了 Mobius 变换，这是 Gibbons 论文中简要描述的一种线性变换。

【参考方案9】：

在我看来，将 e 计算到所需精度的最有效方法是使用以下表示：

e := lim (n -> inf): (1 + (1/n))^n

特别是如果您选择n = 2^x，您可以只用 x 次乘法计算效力，因为：

a^n = (a^2)^(n/2), if n % 2 = 0

【参考方案10】：

二进制拆分方法非常适合模板元程序，它产生一个类型，该类型表示对应于 e 的近似值的有理数。 13 次迭代似乎是最大值 - 任何更高的迭代都会产生“整数常量溢出”错误。

#include <iostream>
#include <iomanip>

template<int NUMER = 0, int DENOM = 1>
struct Rational

    enum NUMERATOR = NUMER;
    enum DENOMINATOR = DENOM;

    static double value;
;

template<int NUMER, int DENOM>
double Rational<NUMER, DENOM>::value = static_cast<double> (NUMER) / DENOM;

template<int ITERS, class APPROX = Rational<2, 1>, int I = 2>
struct CalcE

    typedef Rational<APPROX::NUMERATOR * I + 1, APPROX::DENOMINATOR * I> NewApprox;
    typedef typename CalcE<ITERS, NewApprox, I + 1>::Result Result;
;

template<int ITERS, class APPROX>
struct CalcE<ITERS, APPROX, ITERS>

    typedef APPROX Result;
;

int test (int argc, char* argv[])

    std::cout << std::setprecision (9);

    // ExpType is the type containing our approximation to e.
    typedef CalcE<13>::Result ExpType;

    // Call result() to produce the double value.
    std::cout << "e ~ " << ExpType::value << std::endl;

    return 0;

另一个（非元程序）模板化变体将在编译时计算一个双近似值 e。这个没有迭代次数的限制。

#include <iostream>
#include <iomanip>

template<int ITERS, long long NUMERATOR = 2, long long DENOMINATOR = 1, int I = 2>
struct CalcE

    static double result ()
    
        return CalcE<ITERS, NUMERATOR * I + 1, DENOMINATOR * I, I + 1>::result ();
    
;

template<int ITERS, long long NUMERATOR, long long DENOMINATOR>
struct CalcE<ITERS, NUMERATOR, DENOMINATOR, ITERS>

    static double result ()
    
        return (double)NUMERATOR / DENOMINATOR;
    
;

int main (int argc, char* argv[])

    std::cout << std::setprecision (16);

    std::cout << "e ~ " <<  CalcE<16>::result () << std::endl;

    return 0;

在优化的构建中，表达式 CalcE<16>::result () 将被实际的双精度值替换。

可以说两者都非常有效，因为它们在编译时计算 e :-)

【讨论】：

OVERKILL...不需要这么复杂

【参考方案11】：

@nico 回复：

...比级数的泰勒展开式计算“更快”，即：

e = 1/0！ + 1/1！ + 1/2！ + ...

或

1/e = 1/0！ - 1/1！ + 1/2！ - 1/3！ + ...

以下是在代数上提高牛顿法收敛性的方法：

https://www.researchgate.net/publication/52005980_Improving_the_Convergence_of_Newton's_Series_Approximation_for_e

它们是否可以与二进制拆分结合使用以加快计算速度，这似乎是一个悬而未决的问题。尽管如此，这里有一个来自 Damian Conway 使用 Perl 的示例，它说明了这种新方法在直接计算效率方面的改进。在标题为“? 用于估计”的部分中：

http://blogs.perl.org/users/damian_conway/2019/09/to-compute-a-constant-of-calculusa-treatise-on-multiple-ways.html

（此评论太长，无法在 2010 年 6 月 12 日 10:28 回复回复）

【参考方案12】：

来自wikipedia 将 x 替换为 1

【讨论】：

“由于许多除法运算，将常数 e 表示为无限级数之和的标准表示对于计算非常低效。”

@duffymo：什么是更好的解决方案？

没有说我有更好的解决方案。我正在回复 tur1ng 的评论。