SICP по-русски

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)