Haskell大数计算

Question

我正在尝试用大数字进行一些计算

?: let r  = 291381631919914084
?: let t = 1165526527679656343
?: sqrt(4 * r * r - 4 * r + 1 + 8 * t) - 2 * r + 1
1.0

答案应该8.0000... 是我应该用于此类计算的包吗？还是我应该在前奏中做些什么？

Answer 1

正确的答案确实非常接近8.0.您遇到了数值精度问题:使用IEEE 754("双精度")二进制64格式计算平方根,其53位精度不足以在此给出准确的结果.

更详细:真正的价值sqrt(4 * r * r - 4 * r + 1 + 8 * t)是50个有效数字:

582763263839828175.00000000000000000686385063746811

与该数量最接近的可表示的IEEE 754 binary64值是:

582763263839828224.0

......与真实价值相差约49.0.同样,该值2*r在转换为浮点时会失去精度.

您可能想通过提高精度来解决这个问题,但是在数值工作中经常会发生这种情况,在这种情况下,最好重新编写算法以避免(或至少改善)数值问题.您计算的值是形式sqrt(a * a + b) - a(带a = 2 * r - 1和b = 8 * t).该数量可以在形式被重写b / (sqrt(a * a + b) + a),并(假设两个a和b为正),后者表达将给出更精确的结果.

这是一个快速演示,两个表达式给出相同的结果.

Prelude> let a = 43
Prelude> let b = 7
Prelude> sqrt(a * a + b) - a
8.131845707602992e-2
Prelude> b / (sqrt(a * a + b) + a)
8.131845707603225e-2

我们正在使用的较小值a和b,所以数值的问题是没有那么糟糕,但请注意,仍然在最后4位数字的差异.(这里的确切值是0.08131845707603225000568393232263645035位有效数字.)

并使用这个表达式与您的值:

Prelude> let r = 291381631919914084
Prelude> let t = 1165526527679656343
Prelude> let a = 2*r - 1; b = 8*t in b / (sqrt(a*a+b) + a)
8.0

至于其他的应答者所指出的那样,答案是不准确 8.0的,但8.0 就是最接近的IEEE 754 binary64浮点值,真正的答案.

Answer 2

我相信8也不是正确的答案; 你给出的数字不是正方形:

Math.NumberTheory.Powers.Squares> r  = 291381631919914084
Math.NumberTheory.Powers.Squares> t = 1165526527679656343
Math.NumberTheory.Powers.Squares> isSquare (4*r*r - 4*r + 1 + 8*t)
False

但是,如果你想要的话,你可以得到这个答案:

Math.NumberTheory.Powers.Squares> integerSquareRoot (4*r*r - 4*r + 1 + 8*t) - 2*r + 1
8

所述arithmoi包提供这些功能.

或者,您可以根据需要获得确切答案的数字:

Data.Number.CReal> sqrt (4*r*r - 4*r + 1 + 8*t) - 2*r + 1 :: CReal
8.0000000000000000068638506374681082902485

该数字包提供了这种类型.