Решение упражнения 1.45 из SICP
В первую очередь запишем процедуру, вычисляющую корень n-ой степени из заданного числа, используя заданное количество торможений усреднением:
(define (nth-root n x damp-count) (define (f y) (/ x (power y (- n 1)))) (fixed-point ((repeated average-damp damp-count) f) 1.0))
Здесь используются следующие ранее упоминавшиеся вспомогательные процедуры:
(define (power a b) (exp (* b (log a))))
(define tolerance 0.00001)
(define (fixed-point f first-guess) (define (close-enough? v1 v2) (< (abs (- v1 v2)) tolerance)) (define (try guess) (let ((next (f guess))) (if (close-enough? guess next) next (try next)))) (try first-guess))
(define (average a b) (/ (+ a b) 2.0))
(define (average-damp f) (lambda (x) (average x (f x))))
(define (repeated f n) (if (> n 1) (compose (repeated f (- n 1)) f) f))
Будем проводить эксперименты с корнями из двойки. Сначала установим количество торможений усреднением равным 1:
> (nth-root 2 2 1) 1.4142135623746899 > (nth-root 3 2 1) 1.259923236422975 > (nth-root 4 2 1) ...
На корне 4-ой степени процесс зависает. Используем для 4-ой и последующих степеней два торможения усреднением:
> (nth-root 4 2 2) 1.189207115002721 > (nth-root 5 2 2) 1.1486967244204176 > (nth-root 6 2 2) 1.1224648393618204 > (nth-root 7 2 2) 1.1040857488809646 > (nth-root 8 2 2) ...
Вновь зависаем на 8-ой степени. Начинает вырисовываться закономерность, что увеличивать количество торможений усреднением нужно на степенях двойки, то есть гипотеза такова: для сходимости нашего процесса вычисления корня n-ой степени нужно использовать [log2n] торможений усреднением. Проверяем:
> (nth-root 8 2 3) 1.090507732665258 > (nth-root 15 2 3) 1.0472983541977883 > (nth-root 16 2 3) ... > (nth-root 16 2 4) 1.0442737824274144
Похоже, наша гипотеза подтверждается.
Таким образом окончательный вариант, выведенный, впрочем, чисто эмпирически, выглядит таким образом:
(define (nth-root-empirical n x) (define (f y) (/ x (power y (- n 1)))) (define damp-count (floor (log2 n))) (fixed-point ((repeated average-damp damp-count) f) 1.0))
(define (log2 x) (/ (log x) (log 2)))
Небольшая проверка его работы:
> (nth-root-empirical 4 81) 3.000000000000033 > (nth-root-empirical 2 81) 9.0 > (nth-root-empirical 10 (power 2 10)) 2.0000011830103324 > (nth-root-empirical 100 (power 2 100)) 2.0000032790812043
Что ж, все в порядке.
Comments
Pingback from SICP по-русски » Blog Archive » Сколько программистов нужно, чтобы возвести число A в степень B
Date: January 12, 2008, 11:39 pm
[…] решении упражнения 1.45 я записал процедуру power, возводящую вещественное […]
Comment from thror
Date: February 8, 2008, 9:45 pm
Ну… вообще-то есть достаточно строгое математическое рассуждение, позволяющее получить вашу гипотезу: “для сходимости процесса вычисления корня n-ой степени нужно использовать [log2(n)] торможений усреднением”.
Постараюсь привести здесь его начало, позволяющее выдвинуть такую гипотезу (есть удивительное доказательство ее, но оно слишком громоздко, чтоб уместиться здесь =)).
Итак:
Обозначим корень n-ой степени из a так: x = sqrt_n(a), и a в степени n вот так: x = a**n. Логарифм по основанию 2 от n так log2(n).
Используя метод Ньютона в случае x = sqrt_n(a) имеем
f(x) = x**n – a
функция, нуль которй и есть искомый корень n-ой степени из a. Преобразование Ньютона для нее имеет очевидно вид:
F(x) = (1 – 1/n)x + (1/n)a/(x**(n-1)) ——– [1]
[особое внимание следует обратить на коэффициенты (1-1/n) и (1/n)]
Рассмотрим теперь многократное торможение усреднением для итеративного процесса:
x -> a/(x**(n-1))
k-кратное торможение усреднением дает следующий итеративный процесс:
x -> (1-1/w)x+(1/w)a/(x**(n-1)), w=2**k ——– [2]
[особое внимание следует обратить на коэффициенты (1-1/w) и (1/w)]
Сравните теперь формулы [1] и [2].
Вот теперь, надеюсь, смысл гипотезы на поверхности, а именно, для сходимости итеративного процесса [2] он должен удовлетворять соотношению:
1/w = n
или, что эквивалентно,
k >= ceil(log2(n))
[функия ceil(x) определяется стандартно, например см. Кнут, т.1, Основные Алгоритмы, п. 1.2.4]
Очевидно, наименьшее подходящее k=ceil(log2(n)).
Впрочем, эти рассуждения хоть и подводят нас к данной гипотезе, выводя ее на поверхность, но при этом ни слова не говорят о том, почему все же она верна.
Думаю это будет хорошим упражнением для многих читателей =)
–thror
Comment from thror
Date: February 8, 2008, 9:58 pm
Ключом к доказательству является тот факт, что если модифицировать метод Ньютона, внеся определенные изменения в преобразование Ньютона (какое–не скажу =)) то его сходимость не изменится.
При этом k-кратное торможение усреднением можно представить именно таким модифицированным методом Ньютона.
Ну, вроде многие это даже на уроках проходили =)
Comment from Sergey Khenkin
Date: February 9, 2008, 12:25 am
Хм, очень даже любопытно.
Не совсем понятна мысль:
“для сходимости итеративного процесса [2] он должен удовлетворять соотношению: 1/w = n”.
Comment from thror
Date: February 9, 2008, 10:56 am
А, да, фигня получилась =)
1 / w меньше или равно 1 / n
Comment from thror
Date: February 9, 2008, 2:29 pm
Или w больше или равно n, а т.к. w=2**k, то k больше или равно ceil(log2(n))…
На самом деле рассмотрим преобразование Ньютона для функции f(x):
F(x) = x – f (x) / f ‘ (x)
Модифицируем преобразование таким способом:
G(x) = x – f(x) / (q * f ‘ (x))
Тогда для f(x) = x**n – a получим модифицированное преобразование такое:
G(x) = F(x) = (1 – 1/qn)x + (1/qn)a/(x**(n-1)) ——– [3]
Сравним полученную формулу [3] с формулой [2], и выберем множитель q из условия:
w = qn
Следовательно:
q = w / n
или
q = 2**k / n
Так вот, сходимость модифицированного метода Ньютона не нарушится, когда q больше или равно 1, или что тоже самое, w больше или равно n, или k больше или равно ceil(log2(n)).
Вот собственно и весь секрет этого упражнения.
–thror
Comment from gorilych
Date: September 29, 2008, 2:52 pm
Введём более общую функцию усреднения которая переводит f(x) в 1/m[f(x) + (m-1)x]:
(define (average-damp-m f m) (lambda (x) (/ (+ (f x) (* (- m 1) x)) m)))
Очевидно, что простая функция усреднения average-damp аналогична average-damp-m с параметром m = 2, а (repeated average-damp damp-count) это то же самое, что и average-damp-m при m = 2^(damp-count)
Наш метод “улучшения” создаёт (должен создавать) последовательность чисел, приближающихся к корню n-ной степени из x.
Каждое следующее приближение равно функции улучшения, взятой от предыдущего приближения. В качестве функции улучшения у нас выступает 1/m[x/y^(n-1) + (m-1)y].
Легко показать, что если отклонение предыдущего приближения от решения x^(1/n) составляет малую величину эпсилон (eps/x^(1/n) много меньше 1), то отклонение улучшенного приближения будет приблизительно равно eps*|m-n|/m.
Соответственно, при |m-n|/m больше 1 ряд расходится, при |m-n|/m меньше 1 ряд сходится, при равенстве 1 получаются бесконечные осцилляции вокруг решения. Уточнённое условие сходимости – m > n/2, наилучшая сходимость достигается при m=n.
Comment from gorilych
Date: September 29, 2008, 3:14 pm
Используя новую функцию усреднения достаточно просто создать эффективную функцию вычисления корня n-ной степени:
(define (nth-root n x) (fixed-point (average-damp-m (lambda (y) (/ x (expt y (- n 1)))) n) 1.0))
Comment from potters-owl
Date: July 6, 2010, 3:33 pm
Любопытно, как ув. thror получил
x -> (1-1/w)x+(1/w)a/(x**(n-1)), w=2**k
для многократного торможения усреднением. У меня получается что угодно, кроме этого 🙂
Жаль, что не указано средств обратной связи.
Comment from potters-owl
Date: July 6, 2010, 6:13 pm
впрочем, оказалось достаточно не подставлять сразу у -> x/ y ** (n-1) в формулу торможения усреднением и не пытаться вывести формулу для k-ого торможения, “протаскивая” попутно это отображение n-ого корня, а рассмотреть многократное применение торможения “самого по себе”, как формула сразу стала видна. ЕЩе один довод в пользу подхода, используемого в этом и предыдущем упражнениях )
Comment from Sergk
Date: December 12, 2014, 2:05 pm
Подскажите пожалуйста, первым что пришло в голову написал:
(define (nth-root x n)
(fixed-point (lambda (y) (average y (/ x (accumulate * y identity 1 inc n)))) 1.0))
и тут же не зависимо от степени и без репитов, функция сходится, но потом через часы я всё переписал и нашёл причину + сейчас у вас подсмотрел (всё дело в power), дак вот в чём вопрос: я правильно понимаю, функция сходится из-за того, что результат accumulate это константа,а не функция которая должна получиться?
Comment from Irv
Date: January 14, 2016, 11:58 am
для степеней, не являющимися степенями двойки можно делать на одно усреднение меньше, за счет использования ceiling вместо floor (округление в меньшую сторону, вместо большей) при округлении log2(n).
(define (compose f g) (lambda(x) (f (g x)))) (define (repeat f times) (define (iter i g) (if (= i times) g (iter (+ 1 i) (compose f g)))) (iter 1 f)) (define (average-dump f) (lambda(x) (/ (+ (f x) x) 2))) (define (fixed-point f guess) (define (try x) (let [(y (f x)) (d 0.0000001)] (define (good? x) (> d (abs (- x y)))) (if (good? x) y (try y)))) (try guess)) (define (fixed-point-of-transform f t guess) (fixed-point (t f) guess)) (define (power a b) (exp (* b (log a)))) (define (integer-log k n) (ceiling (/ (log n) (log k)))) (define (root x n) (fixed-point-of-transform (lambda(y) (/ x (power y (- n 1)))) (repeat average-dump (integer-log 2 n)) 1.0)) (root 9 2) (root 27 3) (root 81 4) (root 6561 8) (root 43046721 16) (root 1853020188851841 32)
Write a comment