Квадратные корни по методу Ньютона

Question

Следующая программа Scheme реализует метод Ньютона для вычисления квадрата root числа:

(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
[ Зависание навсегда… ]

Как мы видим, эта наивная программа не работает хорошо:

• для small чисел (ниже x = 1e-2), tolerance 0,001 слишком велико;
• для больших чисел (от x = 1e13), программа зависает навсегда.

Обе проблемы могут быть решены путем переопределения процедуры good-enough? следующим образом:

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

Но это решение не является целью моего поста. Вместо этого я хотел бы понять, почему наивная программа терпит неудачу так, как она терпит неудачу для больших чисел.

Я не читал Стандарт IEEE для арифметики с плавающей точкой c (IEEE 754 ), но, насколько я понимаю, числа с плавающей точкой не могут представлять все вещественные числа и имеют абсолютную точность , которая очень высока для небольших чисел и очень низка для больших чисел (кажется, это цифра из Википедии чтобы подтвердить это). Другими словами, маленькие числа с плавающей точкой являются плотными, а большие числа с плавающей точкой редкими. Следствием этого является то, что наивная программа будет зависать вечно, захваченная бесконечной рекурсией, если guess еще не достиг диапазона допуска и если (improve guess) не может улучшить guess больше, потому что расстояние между новыми guess, а старый guess ниже абсолютной точности старого guess (поэтому новый guess - это то же самое число с плавающей точкой, что и старый guess, что означает, что фиксированная точка (improve guess) имеет было достигнуто).

Чтобы гарантировать, что для данного x наивная программа вернется, мне кажется, что этот предикат должен содержать:

допуск> absolute_precision (sqrt (x)).

Если мы выберем tolerance из 0,001 = 1e-3, это означает, что абсолютная точность sqrt (x) должна быть меньше 1e-3 , Следовательно, согласно приведенному выше рисунку из Википедии для двоичных чисел с плавающей запятой 64, sqrt (x) должно быть меньше 1e13, и поэтому x должно быть меньше 1e26.

Вот мои вопросы:

1. Верен ли мой предикат?
2. Почему программа навсегда зависает от x = 1e13 вместо ожидаемого x = 1e26?
3. Почему программа висеть вечно с x в {1e13, 1e15, 1e17,…, 1e49} и любых x больше, чем 1e49, но все еще возвращаться с x в {1e14, 1e16, 1e18,…, 1e48}?

Примечание. - Я использую интерпретатор Chibi Scheme на 64-битном MacBook Pro, поэтому формат IEEE 754 binary64.

Answer 1

Когда вы возводите в квадрат 3162277.6601683795, вы получаете 10000000000000.002 (в двоичном виде это бесконечно, поэтому теряет точность по сравнению с аргументом). Если вы конвертируете 10000000000000002 в гекс, вы получаете 2386F26FC10002, что составляет 4x13 клев, а 2 использует 2 бита. В float значение msb всегда равно 1 и опущено, поэтому я предполагаю, что используется 53 бита и пропускается много десятичных знаков, поскольку именно так работает плавающая точка. Вы не можете приблизиться к 10000000000000 по допуску 0.001 с любой стороны

Я думаю, вы могли бы сравнить предположение с предыдущим предположением и посмотреть, не превышает ли оно допуск. например.

(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.

Теперь, если вы посмотрите на предположения, которые он делает, он будет прыгать довольно сильно, делая новое предположение всегда намного лучше, поэтому только когда предположение близко к фактическому результату, оно будет близко к предыдущему предположению, и я предполагаю, что вы получили бы такую ​​же хорошую оценку, если бы ваша терпимость была просто возведена в квадрат от предыдущей

Я согласен с тем, что фиксированная толерантность не годится, и у нас должно было быть что-то, что смотрело бы на догадки, разделенные между собой. Когда они близки, результат всегда будет близок к 1, а в то время как он будет выключен, он будет близок к 2 или 1/2 в зависимости Простое сравнение без допусков может привести к бесконечному значению l oop, когда improve будет перепрыгивать через отметку и никогда не совпадать, поэтому допуск был установлен там по причине.