英文:
how to fix the precision for calculating napier's constant
问题
我尝试使用无限数组计算自然常数 e。
a = 1.0:0.1:[y/10 | (x,y) <- zip a (tail a)]
b = map (\n -> (1+n)**(1/n)) a
前十几个数的结果是正常的，但对于较小的数不正常：
Prelude> b!!2
2.7048138294215285
Prelude> b!!3
2.7169239322355936
Prelude> b!!4
2.7181459268249255
Prelude> b!!5
2.718268237192297
Prelude> b!!6
2.718280469095753
Prelude> b!!500
1.0
有什么办法可以修复吗？
英文:
I tried to calculate the napier's constant e using the infinite array.
a = 1.0:0.1:[y/10 | (x,y) <- zip a (tail a)]
b = map (\n -> (1+n)**(1/n)) a
The result of first dozen of the number was ok but not for smaller numbers;
Prelude> b!!2
2.7048138294215285
Prelude> b!!3
2.7169239322355936
Prelude> b!!4
2.7181459268249255
Prelude> b!!5
2.718268237192297
Prelude> b!!6
2.718280469095753
Prelude> b!!500
1.0
Any idea how to fix it?
答案1
得分: 2
您的方法等同于更容易理解的方法：
b = [ (1+1/n)^n | i<-[0..], let n=10^i ]
请注意，每个项都是单独计算的，因此没有必要将它们计算为无限列表 - 您可以直接计算具有足够精度的那个项，并跳过其他所有项。
您观察到的问题是，这个算法在数值上不稳定；在实践中不应该使用。正如chi所指出的问题，浮点运算中的 1+n 最终会评估为精确的一，当 n 足够小。并且将1升到任何幂次方总是得到1。
您可能认为可以通过切换到精确算术而不是浮点算术来解决此问题。事实上，您可以理论上使用类型 Rational 来评估上面的表达式，而这将正确收敛。但是您会发现这是不可行的，因为分数随着它们被升到某个幂次方而变得越来越大，计算在您获得任何优势之前就停滞了。
此外，这只能工作是因为我将它从 n=10^^(-i) 然后到 **(1/n) 改为了整数 n=10^i，这允许使用 ^ 运算符。** 支持分数参数，但对于有理数来说是不可能的（分数指数对应于根，通常是无理数）。
真正讽刺的是 ** 是基于 指数 实现的：
x**y = exp(log x*y)
因此，您的原始版本实际上会评估如下：
exp(log(1 + n)/n)
– 再次不稳定（您会得到 log(1+n) == 0），但此时只需在解析上取极限就可以了，即
exp((0 + n + 𝑛²)/n)
→ exp(n/n)
≡ exp 1
– 嗯，当然这是欧拉/纳皮尔常数... 但它相当自我引用，不是吗？
如果您想以实际有意义的方式计算该常数，您应该使用级数展开而不是这种方法。
英文:
Your approach is equivalent to the much more understandable
b = [ (1+1/n)^n | i<-[0..], let n=10^i ]
Note that each term is computed individually, so there's not much point in computing them as an infinite list – you might as well just compute the one that has sufficient accuracy right away and skip all the others.
The problem you're observing is that this algorithm is numerically unstable; it should not be used in practice. The problem, as chi remarked, is that 1+n in floating-point arithmetic eventually evaluates to exactly one, with small enough n. And raising 1 to any power gives always 1.
You might think this could be circumvented by switching to exact arithmetic instead of floating-point. Indeed you can theoretically evaluate the above with the type Rational, and that would properly converge, however you'll find that it's utterly infeasible because the fractions get bigger and bigger as they're raised to a power, and the computation grinds to a halt long before you gain any advantage over the float version.
Also, this only works because I changed it from n=10^^(-i) and then **(1/n) to the integral n=10^i, which allows using the ^ operator. ** supports fractional arguments, but that's not possible for rational numbers (fractional exponents correspond to roots, which are generally irrational).
The really ironic thing is that ** is implemented in terms of exponentials:
x**y = exp(log x*y)
So your original version actually evaluates like
exp(log(1 + n)/n)
– again, unstable (you get log(1+n) == 0), but at this point it would have been trivial to just take the limit analytically, which is
exp((0 + n + 𝓞(n²))/n)
→ exp(n/n)
≡ exp 1
– and, well, sure that is the Euler/Napier constant... but it's rather self-referential, isn't it?
If you want to compute that constant in an actually meaningful way, you should use a series expansion instead.
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论