牛顿法求平方根 [英] Square roots by Newton’s method

问题描述

以下Scheme程序实现了牛顿方法，用于计算数字的平方根:

The following Scheme program implements Newton’s method for computing the square root of a number:

(import (scheme small))

(define (sqrt x)
  (define (sqrt-iter guess)
    (if (good-enough? guess)
      guess
      (sqrt-iter (improve guess))))
  (define (good-enough? guess)
    (define tolerance 0.001)
    (< (abs (- (square guess) x)) tolerance))
  (define (improve guess)
    (if (= guess 0) guess (average guess (/ x guess))))
  (define (average x y)
    (/ (+ x y) 2))
  (define initial-guess 1.0)
  (sqrt-iter initial-guess))

(display (sqrt 0)) (newline)
(display (sqrt 1e-12)) (newline)
(display (sqrt 1e-10)) (newline)
(display (sqrt 1e-8)) (newline)
(display (sqrt 1e-6)) (newline)
(display (sqrt 1e-4)) (newline)
(display (sqrt 1e-2)) (newline)
(display (sqrt 1e0)) (newline)
(display (sqrt 1e2)) (newline)
(display (sqrt 1e4)) (newline)
(display (sqrt 1e6)) (newline)
(display (sqrt 1e8)) (newline)
(display (sqrt 1e10)) (newline)
(display (sqrt 1e12)) (newline)
(display (sqrt 1e13))

输出:

0.03125
0.031250000010656254
0.03125000106562499
0.03125010656242753
0.031260655525445276
0.03230844833048122
0.10032578510960605
1.0
10.000000000139897
100.00000025490743
1000.0000000000118
10000.0
100000.0
1000000.0
[永远挂…]

0.03125
0.031250000010656254
0.03125000106562499
0.03125010656242753
0.031260655525445276
0.03230844833048122
0.10032578510960605
1.0
10.000000000139897
100.00000025490743
1000.0000000000118
10000.0
100000.0
1000000.0
[Hanging forever…]

如我们所见，这个幼稚的程序表现不佳:

As we can see, this naive program does not perform well:

• 对于小数字(在x = 1e-2以下)，tolerance 0.001太大；
• 对于大个数字(从x = 1e13开始)，该程序将永远挂起.

• for small numbers (below x = 1e-2), the tolerance 0.001 is too large;
• for large numbers (from x = 1e13), the program hangs forever.

可以通过重新定义good-enough?过程来解决这两个问题:

Both problems can be solved by redefining the good-enough? procedure like this:

(define (good-enough? guess)
  (= (improve guess) guess))

但是这种解决方案不是我发帖的目的.相反，我想了解为什么朴素的程序会以失败的方式失败，而不是大量失败.

But this solution is not the purpose of my post. Instead, I would like to understand why the naive program fails the way it fails for large numbers.

我尚未阅读 IEEE浮动标准点算术(IEEE 754)，但据我了解，浮点数不能表示所有实数，并且具有绝对精度，对于小数而言非常高，对于大数而言非常低(此维基百科图似乎证实了这一点).换句话说，小的浮点数是密集的而大的浮点数是稀疏的.这样的结果是，如果guess尚未达到公差范围，并且如果(improve guess)不能再改善guess，因为新guess并且旧的guess低于旧的guess的绝对精度(因此新的guess是与旧的guess相同的浮点数，这意味着(improve guess)的固定点具有已达到).

I have not read the IEEE Standard for Floating-Point Arithmetic (IEEE 754), but to my understanding, floating-point numbers cannot represent all reals and have an absolute precision that is very high for small numbers and very low for large numbers (this Wikipedia figure seems to confirm this). In other words, small floating-point numbers are dense and large floating-point numbers are sparse. The consequence of this, is that the naive program will hang forever, trapped in an infinite recursion, if guess has not reached the tolerance range yet and if (improve guess) cannot improve the guess anymore because the distance between the new guess and the old guess is below the absolute precision of the old guess (so the new guess is the same floating-point number as the old guess, meaning that a fixed point of (improve guess) has been reached).

为了保证对于给定的x，朴素的程序返回，在我看来，该谓词必须成立:

To guarantee that, for a given x, the naive program returns, it seems to me that this predicate must hold:

tolerance> absolute_precision(sqrt(x)).

如果我们选择0.001 = 1e-3的tolerance，这意味着sqrt(x)的绝对精度应小于1e-3.因此，根据上面的Wikipedia图中的binary64浮点数，sqrt(x)应该小于1e13，因此x应该小于1e26.

If we choose a tolerance of 0.001 = 1e-3, that means that the absolute precision of sqrt(x) should be less than 1e-3. Consequently, according to the Wikipedia figure above for binary64 floating-point numbers, sqrt(x) should be less than 1e13 and therefore x should be less than 1e26.

这是我的问题:

1. 我的谓词正确吗?
2. 为什么程序会从x = 1e13永久挂起，而不是预期的x = 1e26?
3. 为什么程序会在{1e13，1e15，1e17，…，1e49}中和x以及大于1e49的任何x一起永久挂起，而在{1e14，1e16，1e18，…，中仍然以x返回， 1e48}?

1. Is my predicate correct?
2. Why does the program hang forever from x = 1e13 instead of the expected x = 1e26?
3. Why does the program hang forever with x in {1e13, 1e15, 1e17, …, 1e49} and any x greater than 1e49 but still return with x in {1e14, 1e16, 1e18, …, 1e48}?

注意.-我正在使用 Chibi Scheme 解释器在64位MacBook Pro上使用IEEE 754 binary64格式.

Note. — I am using the Chibi Scheme interpreter on a 64-bit MacBook Pro so the IEEE 754 binary64 format.

推荐答案

对3162277.6601683795求平方会得到10000000000000.002(二进制形式是无限的，因此与参数相比，精度降低).如果将10000000000000002转换为十六进制，则会得到2386F26FC10002，它是4x13个半字节，而2使用2位.在float中，msb始终为1，因此被省略，因此我猜测使用了53位，并且跳过了许多小数位，因为这就是浮点的工作方式.您无法从任一侧通过公差0.001来更接近10000000000000

When you square 3162277.6601683795 you get 10000000000000.002 (in binary this is infinite so looses precision compared to the argument). If you convert 10000000000000002 to hex you get 2386F26FC10002 which is 4x13 nibbles and 2 uses 2 bits. In float the msb is always 1 and is omitted so I'm guessing 53 bits is used and lots of decimals are skipped since that is how floating point works. You cannot get closer to 10000000000000 by tolerance 0.001 from either side

我在想，您可以将猜测与先前的猜测进行比较，看看它是否在公差范围内.例如.

I'm thinking you could compare the guess with the previous guess and see if that is below a tolerance. eg.

(define (sqrt x)
  (define (sqrt-iter guess prev-guess)
    (if (good-enough? guess prev-guess)
      guess
      (sqrt-iter (improve guess) guess)))
  (define (good-enough? guess prev-guess)
    (define tolerance 1e-20)
    (< (abs (- guess prev-guess)) tolerance))
  (define (improve guess)
    (if (= guess 0) guess (average guess (/ x guess))))
  (define (average x y)
    (/ (+ x y) 2))
  (define initial-guess 1.0)
  (sqrt-iter initial-guess 0))

这将导致3162277.6601683795.与您的版本不同，当您减去猜测时，丢失的位会使差值变小，因此即使容差为1e-20，它也会为1e270给出答案.

This results in 3162277.6601683795. Unlike your version when you substract the guesses missing bits will make the difference smaller so even with a tolerance of 1e-20 it will produce an answer for 1e270.

现在，如果您查看这些猜测，它将使跳跃幅度变大，从而使新的猜测始终更好，因此只有当该猜测接近实际结果时，它才会接近先前的猜测，我猜您会如果您的容忍度只是前一个的平方，则可以得到很好的估计.

Now if you look at the guesses it makes it will jump quite large making the new guess always much better so only when the guess is close to the actual result it will be close to the previous guess and I'm guessing you would get as good estimate if your tolerance were just squared from the previous.

我同意固定容忍度不好，我们应该有一些看待彼此猜测的猜测.当它们接近时，结果将始终接近1，而当偏离时，结果将接近2或1/2，具体取决于.仅当没有公差的情况下进行比较时，当improve会在标记上方和下方跳跃并且永远不匹配时，您可能会陷入无限循环，因此将公差放入其中是有原因的.

I agree that the fixed tolerance is no good and we should have had something that looks at the guesses divided with each other instead. When they are close the result will always be close to 1 and while it is off by much it will be close to 2 or 1/2 depending. Just comparing without tolerance might put you in a an infinite loop when the improve will jump over and under the mark and never match so the tolerance was put in there for a reason.

这篇关于牛顿法求平方根的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！