John Harrison - Введение в функциональное программирование (1108517), страница 31
Текст из файла (страница 31)
Точная арифметика вещественных чиселГлава 9. Примеры1 11+ |fx (n + 2) − 2n+2 x| + |fy (n + 2) − 2n+2 y|2 4411 1+ 1+ 1<2 44= 1≤Программная реализация выглядит так:l e t real_add f g n =( f ( n + 2 ) +/ g ( n + 2 ) ) ndiv ( I n t 4 ) ; ;Аналогичные рассуждения могут быть использованы для определения операциивычитания, но проще всего построить её на основе уже введенных функций:#let real_sub f g = real_add f (real_neg g);;real_sub : (num -> num) -> (num -> num) -> num -> num = <fun>Реализация умножения, вычисления обратных чисел и деления потребует в общемслучае несколько больших усилий.
Однако, частные случаи умножения и деленияна целое число существенно легче и при этом достаточно часты. По соображениямэффективности они заслуживают отдельного рассмотрения. Пустьfmx (n) = (mfx (n + p + 1)) ndiv 2p+1 ,где p выбирается так, чтобы 2p ≥ |m|. Доказать корректность этого определениялегко:|fmx (n) − 2n (mx)| ≤=<≤≤1212121212mfx (n + p + 1)− 2n (mx)|2p+1|m||fx (n + p + 1) − 2n+p+1 x|2p+1|m|2p+11 |m|2 2p1= 1.2+|++++Для реализации такого подхода нам потребуется функция вычислениясоответствующего p. Не слишком изощрённое, но вполне подходящее определениеможет выглядеть так:let log2 =l e t rec l o g 2 x y =i f x </ I n t 1 then ye l s e l o g 2 (quo_num x ( I n t 2 ) ) ( y + 1 ) infun x −> l o g 2 ( x −/ I n t 1 ) 0 ; ;С учётом сказанного выше, операция умножения на целое число принимает вид:136Глава 9. Примеры9.3.
Точная арифметика вещественных чиселlet real_intmul m x n =l e t p = l o g 2 ( abs_num m) inl e t p1 = p + 1 in(m ∗/ x ( n + p1 ) ) ndiv ( I n t 2 ∗∗/ I n t p1 ) ; ;Деление на целое число вводится следующим образом:fx/m (n) = fx (n) ndiv mДля упрощения доказательства корректности этого определения будем считать,что случай m = 0 никогда не может возникнуть, а при m = ±1 результатоперации не изменяет погрешности. В остальных случаях из |fx (n) − 2n x| < 11следует, что |fx (n)/m − 2n x/m| < |m|≤ 12 , откуда, в свою очередь, с учётом|fx (n) ndiv m − fx (n)/m| ≤ 12 получаем требуемое. В итоге, программная реализацияделения такова:let real_intdiv m x n =x ( n ) ndiv ( I n t m) ; ;9.3.4Умножение: общий случайОпределить умножение в общем случае труднее, поскольку погрешностьаппроксимации одного из сомножителей умножается на порядок второго.Следовательно, нам потребуется предварительно оценить порядки сомножителей.Поступим следующим образом.
Предположим, что нам требуется вычислитьвыражение x + y до n-го разряда. Для этого выберем r и s такие, что |r − s| ≤ 1и r + s = n + 2. Таким образом, обе эти величины несколько больше, чем половинатребуемой разрядности. Далее вычислим fx (r) и fy (s), после чего определим p и q —соответствующие «двоичные логарифмы», для которых справедливо |fx (r)| ≤ 2p and|fy (s)| ≤ 2q . Если как p, так и q равны нулю, то легко убедиться, что результатомвычислений также будет 0. В противном случае отметим, что либо p > 0, либо q > 0 —этот факт нам понадобится в дальнейшем.Пустьk = n+q−s+3=q+r+1l = n+p−r+3=p+s+1m = (k + l) − n = p + q + 4.Покажем, что погрешность выражения fxy (n) = (fx (k)fy (l)) ndiv 2m удовлетворяеткритерию аппроксимации, т. е. |fxy (n) − 2n (xy)| < 1.
Если ввести обозначения2k x = fx (k) + δ,2l y = fy (l) + ,где |δ| < 1 и || < 1, то мы получим:|fxy (n) − 2n (xy)| ≤fx (k)fy (l)1+|− 2n (xy)|22m1379.3. Точная арифметика вещественных чисел===≤≤<121212121212Глава 9. Примеры+ 2−m |fx (k)fy (l) − 2k+l xy|+ 2−m |fx (k)fy (l) − (fx (k) + δ)(fy (l) + )|+ 2−m |δfy (l) + fx (k) + δ|+ 2−m (|δfy (l)| + |fx (k)| + |δ|)+ 2−m (|fy (l)| + |fx (k)| + |δ|)+ 2−m (|fy (l)| + |fx (k)| + 1)Отсюда имеем |fx (r)| ≤ 2p , так что |2r x| < 2p + 1. Следовательно, |2k x| < 2q+1 (2p + 1),откуда |fx (k)| < 2q+1 (2p + 1) + 1, т.
е. |fx (k)| ≤ 2q+1 (2p + 1). Аналогично доказываетсясправедливость |fy (l)| ≤ 2p+1 (2q + 1). Таким образом,|fy (l)| + |fx (k)| + 1 ≤ 2p+1 (2q + 1) + 2q+1 (2p + 1) + 1= 2p+q+1 + 2p+1 + 2p+q+1 + 2q+1 + 1= 2p+q+2 + 2p+1 + 2q+1 + 1.Дляполучениятребуемоймаксимальнойпогрешностивведёмm−1ограничение |fy (l)| + |fx (k)| + 1 ≤ 2, либо, разделив на 2 и учитывая дискретностьмножества целых чисел,2p+q+1 + 2p + 2q < 2p+q+2 .В свою очередь, данное отношение может быть записано как (2p+q + 2p ) + (2p+q +2q ) < 2p+q+1 + 2p+q+1 .
Его справедливость следует из упомянутого ранее факта, чтолибо p > 0, либо q > 0. Таким образом, мы обосновали следующее определение:l e t real_mul x y n =l e t n2 = n + 2 inl e t r = n2 / 2 inl e t s = n2 − r inl e t xr = x ( r )and ys = y ( s ) inl e t p = l o g 2 xrand q = l o g 2 ys ini f p = 0 & q = 0 then I n t 0 e l s elet k = q + r + 1and l = p + s + 1and m = p + q + 4 in( x ( k ) ∗/ y ( l ) ) ndiv ( I n t 2 ∗∗/ I n t m) ; ;9.3.5Обратные числаНашим следующим шагом будет реализация вычисления обратных чисел.
Чтобыполучить любую верхнюю оценку обратного числа, не говоря уж о хорошем егоприближении, потребуется оценить аргумент снизу. В общем случае для этого не138Глава 9. Примеры9.3. Точная арифметика вещественных чиселсуществует лучшего способа, чем вычисление аргумента с возрастающей точностью,пока we can bound it away from zero. Следующая лемма служит обоснованием этойпроцедуры.Лемма 9.1 Пусть 2e ≥ n + k + 1, |fx (k)| ≥ 2e и |fx (k) − 2k x| < 1, где fx (k) — целое,а e, n и k — натуральные числа. Если мы определимfy (n) = 2n+k ndiv fx (k),то получим |fy (n) − 2n x−1 | < 1, т.
е. требуемую верхнюю оценку погрешности.Доказательство: Доказательство этой леммы достаточно утомительно и будетприведено здесь не полностью, а лишь в виде основных соображений. Если |fx (k)| >2e , то результат очевидно следует из простых рассуждений. Округление врезультате выполнения операции ndiv даёт погрешность, не превышающую 21 ,итоговая погрешность меньше 12 . Если |fx (k)| = 2e , но при этом n + k ≥ e,то можно пренебречь тем, что второй компонент погрешности может бытьвдвое большим, т. е. меньшим 1 — в этом случае ошибка округления будетотсутствовать потому, что fx (k) = ±2e делится на 2n+k нацело.
(Воспользуемсяфактом, что 2e − 1 ≤ 2e−1 , поскольку при 2e ≥ n + k + 1 значение e не можетравняться нулю.) Наконец, при |fx (k)| = 2e и n + k < e, мы имеем |fy (n) − 2n x1 | < 1,поскольку |fy (n)| ≤ 1 и 0 < |2n x1 | < 1, а знаки этих величин совпадают. Предположим, что нам требуется найти число, обратное x, с точностью n. Преждевсего, вычислим fx (0).
Нам потребуется рассмотреть два случая:1. Если |fx (0)| > 2r для некоторого натурального числа r, то выберем наименьшеенатуральное число k (которое может быть и равным нулю) такое, что 2r +k ≥ n + 1, и положим e = r + k. Результатом вычислений в этомслучае будет 2n+k ndiv fx (k). Легко убедиться, что условия, требуемые леммой,выполняются.
Поскольку |fx (0)| ≥ 2r +1, мы имеем |x| > 2r , так что |2k x| > 2r+k .Это значит, что |fx (k)| > 2r+k − 1, а отсюда (учитывая целочисленность fx (k))получаем требуемое соотношение |fx (k)| ≥ 2r+k = 2e . Условие 2e ≥ n = k + 1легко проверить. Отметим, что из r ≥ n непосредственно следует корректностьаппроксимации fy (n) = 0.2. При |fx (0)| ≤ 1 воспользуемся функцией ‘msd’, возвращающей наименьшее pтакое, что |fx (p)| > 1.
Отметим, что при x = 0 произойдёт зацикливание.Положив e = n + p + 1 и k = e + p, определим результат операциикак 2n+k ndiv fx (k). В этом случае также справедливы условия, требуемыелеммой. Так как |fx (p)| ≥ 2, мы имеем |2p x| > 1, т. е. |x| > 21p .Следовательно |2k x| > 2k−p = 2e , откуда |fx (k)| > 2e − 1, т. е. |fx (k)| ≥ 2e .Реализацию начнём с функции msd:l e t msd =l e t rec msd n x =i f abs_num ( x ( n ) ) >/ I n t 1 then n e l s e msd ( n + 1 ) x inmsd 0 ; ;1399.3.
Точная арифметика вещественных чиселГлава 9. Примерыпосле чего воплотим изложенные выше теоретические рассуждения в виде простойпрограммы:let real_inv x n =l e t x0 = x ( 0 ) inlet k =i f x0 >/ I n t 1 thenl e t r = l o g 2 x0 − 1 inl e t k0 = n + 1 − 2 ∗ r ini f k0 < 0 then 0 e l s e k0elsel e t p = msd x inn + 2 ∗ p + 1 in( I n t 2 ∗∗/ I n t ( n + k ) ) ndiv ( x ( k ) ) ; ;В итоге определение операции деления становится тривиальным:l e t r e a l _ d i v x y = real_mul x ( r e a l _ i n v y ) ; ;9.3.6Отношения порядкаПримечательным свойством отношений порядка является то, что их вычислениев общем случае алгоритмически неразрешимо. В основе такого вывода лежитневозможность установить в общем случае равенство данного числа нулю. Еслипоследовательное вычисление с увеличением требуемой точности продолжаетвыдавать 0, это ещё не служит гарантией того, что в дальнейшем мы не получимненулевой результат.5 Если значение x не равно нулю, поиск первого ненулевогоразряда числа когда-либо завершится, но в случае x = 0 он будет длиться вечно.Принимая во внимание сказанное выше, несложно реализовать отношенияпорядка.
Для определения взаимного порядка чисел x и y достаточно найти n такое,что |xn − yn | ≥ 2. Например, для xn ≥ yn + 2 мы имеем2n x > xn − 1 ≥ yn + 1 > 2n y.откуда делаем вывод, что x > y. Прежде всего, приведем общую процедурувычисления n, после чего все отношения порядка могут быть выражены с еёиспользованием.
Отметим, что единственным способом реализации рефлексивныхотношений будет положить их тождественными соответствующим нерефлексивнымотношениям!5Доказательство алгоритмической неразрешимости этой задачи сводится к проблемезавершимости: пусть f (n) = 1, если некоторая машина Тьюринга завершает работу не более чем заn итераций, и f (n) = 0 в противном случае.140Глава 9.
Примеры9.3. Точная арифметика вещественных чиселlet separate =l e t rec s e p a r a t e n x y =l e t d = x ( n ) −/ y ( n ) ini f abs_num ( d ) >/ I n t 1 then de l s e s e p a r a t e ( n + 1 ) x y inseparate 0 ; ;letletletlet9.3.7real_gtreal_gereal_ltreal_lexxxxyyyy====s e p a r a t e x y >/ I n t 0 ; ;real_gt x y ; ;s e p a r a t e x y </ I n t 0 ; ;real_lt x y ; ;КэшированиеЧтобы протестировать определённые нами функции, потребуется возможностьвывода некоторого приближённого значения вещественного числа в десятеричнойсистеме счисления. Возможности стандартной библиотеки CAML делают эту задачунесложной. Если нам требуется вывести d десятичных знаков числа, положимточность вычислений n такой, чтобы 2n > 10d , т.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
