О компьютерном округлении вещественных чисел к целым числам

AlgART / Тексты

О компьютерном округлении вещественных чисел к целым числам

Даниэль Алиевский

Выражаю благодарность Константину Кобылкину за ценные замечания и дополнения, в частности за первоначально упущенное мной требование 3 к компьютерному представлению вещественных чисел

Требования к компьютерному представлению вещественных чисел
Лемма о представимости целых и полуцелых
Лемма о монотонности сложения
Теорема округления
Практическое выполнение округления по алгоритму «прибавить 0.5 и взять целую часть»
Теорема вычитания
Учет теоремы округления в Java 7

В современных языках программирования для округления вещественного числа x к ближайшему целому обычно используется операция, которую можно записать как ⌊x+0.5⌋, т.е. наибольшее целое число, не превосходящее x+0.5 (иными словами, целая часть x+0.5). Как правило, для этого существует та или иная библиотечная функция. В частности, в языке Java (одном из лучших с точки зрения точной формальной спецификации вычислений) есть функция StrictMath.round(double x), результат которой, согласно документации (во всех версиях языка, предшествовавших Java 7), равен результату выполнения операторов
(long)Math.floor(x + 0.5),
что буквально соответствует описанному выше правилу.

Однако известно, что компьютерные вычисления выполняются неточно. Действительный результат вычислений нередко отличается от того, что дает строгая математическая формула, и из-за этого выводы, которые были бы истинными в случае математически точных формул, иногда оказываются неверными для реальных программ, что приводит к трудноуловимым ошибкам. К счастью, результат операции округления, определенной описанным выше способом, почти точно соответствует аккуратной математической формуле.

Центральная тема данной работы — теорема округления, специфицирующая это «почти». С помощью этой теоремы можно дать четкие прогнозы поведению программ, выполняющих округление описанным выше способом, т.е. прибавлением 0.5 и взятием целой части. Один из примеров выводов, которые можно сделать относительно использования такого округления, — теорема вычитания, приведенная в конце статьи.

Требования к компьютерному представлению вещественных чисел

Мы будем считать, что компьютерная арифметика реализована в соответствии с определенным стандартом, таким как IEEE 754, и что операции выполняются над числами с плавающей точкой определенного типа, такими как типы IEEE 754 “Double precision” (64 бита, тип double в традиционных языках программирования) или “Single precision” (32 бита, тип float в традиционных языках программирования). Как правило, на современных компьютерах действительно приходится иметь дело с типами double или float, соответствующими IEEE 754; в случае языка программирования Java это требование заложено в спецификацию языка. Однако в действительности сформулированные далее теоремы используют лишь часть требований стандарта IEEE 754 и остаются верными для других реализаций вещественных чисел, соответствующих этим требованиям.

Ниже приведены требования к реализации вещественных чисел, которые мы реально используем. В частности, эти требования удовлетворяются для всех двоичных вещественных форматов стандарта IEEE 754.

Конечное вещественное число x может быть точно представлено на компьютере тогда и только тогда, когда оно равно 0 либо имеет вид ±2^q(1+m/2^P), где P — целая положительная константа, m (мантисса) — произвольное целое число из диапазона 0≤m<2^P, q (порядок) — произвольное целое число из некоторого константного диапазона q_min≤q≤q_max (q_min<0, q_max>0).

Иначе говоря, мы говорим о нормализованных двоичных вещественных числах, таких как привычные double или float. Константа P — это число двоичных разрядов, отведенных под хранение мантиссы (кроме старшего единичного разряда, который в таких форматах не хранится); в случае double это 52, в случае float это 23. Заметим, что на самом деле мы требуем нормализованности лишь для упрощения доказательства — при желании нетрудно распространить приведенное далее доказательство и на случай, когда особо малые значения (близкие к 0) могут быть представлены в ненормализованном общем виде ±2^q−Pm.

Числа указанного вида, которые допускают точное представление на компьютере, мы будем далее кратко называть представимыми.

Очевидным следствием данного определения является независимость представимости числа от знака: оба числа +x и −x всегда одновременно представимы или одновременно непредставимы.

Нижняя граница диапазона порядков удовлетворяет условию: q_min<−P.

Практически во всех реализациях вещественных чисел это так. В частности, это верно для чисел double и float стандарта IEEE 754 (где q_min=−1022 и q_min=−126 соответственно).

В действительности нам нужно лишь следующее простое следствие из этого условия: все числа вида ±m/2^P+1, 0≤m<2^P, являются представимыми. Действительно, при 2^P−1≤m<2^P (диапазон 0.25≤|x|<0.5) такие числа суть ±m/2^P+1=±2⁻²(1+(2m−2^P)/2^P), при 2^P−2≤m<2^P−1 (диапазон 0.125≤|x|<0.25) такие числа суть ±m/2^P+1=±2⁻³(1+(4m−2^P)/2^P) и т.д. Завершает эту последовательность диапазонов единственное значение m=1, при котором 2^P−P≤m<2^P−P+1, и мы имеем ±1/2^P+1=±2^−(P+1)(1+0/2^P) — при q_min<−P это число все еще представимо. Оставшийся не рассмотренным вариант m=0 (т.е. x=0) является исключением и оговорен в требовании 1 отдельно.

Верхняя граница диапазона порядков удовлетворяет условию: q_max>P.

Как и в случае 2-го требования, практически во всех реализациях вещественных чисел это так. В частности, это верно для чисел double и float стандарта IEEE 754 (где q_max=1023 и q_max=127 соответственно).

Это условие позволяет гарантировать представимость достаточно больших по абсолютной величине целых чисел: см. далее лемму о представимости целых и полуцелых.

На самом деле это условие можно было бы без вреда для дела ослабить до q_max≥P, но это привело бы к несколько более громоздким формулировкам последующих утверждений; в частности, при q_max=P числа ±2^P+1 оказались бы непредставимыми.

При компьютерном сложении двух представимых вещественных чисел a и b получается представимое число, максимально близкое к точной математической сумме a+b. (Если таких представимых чисел два, то получается одно из них; какое именно, для нас неважно.)

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

Результат компьютерного сложения представимых двух вещественных чисел a и b мы будем далее обозначать специальным образом: a⊕b. Обычные знаки математических операций (+, −, *, /) с этого момента будут означать результат абсолютно точного выполнения математической операции (который не всегда совпадает с результатом фактических компьютерных вычислений).

Из правила 4 необходимо сделать одно исключение, а именно, случай, когда |a+b| больше максимально представимого вещественного числа, т.е. |a+b|>2^q_max(1+(2^P−1)/2^P). В такой ситуации реальное компьютерное сложение может привести к переполнению, которое детектируется тем или иным способом, например, генерацией исключения или возвратом особого значения, обозначающего бесконечность. В этом специальном случае мы заменим общее требование 4 более частным, а именно: если |a+b|>2^q_max(1+(2^P−1)/2^P) и при этом b=0.5, то компьютерное сложение должно давать a⊕b=a (что вполне естественно в силу требования 3: если a+b по абсолютной величине настолько велико, то значение b=0.5 при сложении пренебрежимо мало и должно быть заменено нулем). Это дополнение на самом деле не очень существенно и нужно лишь для упрощения формулировки теоремы округления (случай |x|>2^P+1).

Очевидное следствие указанного правила: для любых представимых a и b a⊕b=a+b тогда и только тогда, когда математическая сумма a+b является представимой. Действительно, если a⊕b=a+b, то это представимое число: компьютерное сложение, как и все прочие компьютерные операции, по определению возвращает представимые числа. Наоборот, если a+b есть представимое число, то, во-первых, |a+b| никак не может быть больше максимально представимого вещественного числа (ведь a+b представимо, а представимость не зависит от знака), а во-вторых, никакое другое представимое число не может быть ближе к a+b, чем оно само. А это значит, что a⊕b=a+b.

Лемма о представимости целых и полуцелых

Данная лемма является тривиальным следствием требований 1 и 3 к реализации компьютерной арифметики, указанных выше. Однако она имеет большое практическое значение, поэтому есть смысл сформулировать ее отдельно.

Лемма о представимости целых и полуцелых:
Если вещественное число x является целым и при этом |x|≤2^P+1, то такое число x является представимым.
Если вещественное число x является полуцелым, т.е. имеет вид k+0.5, где k целое, и при этом |x|≤2^P, то такое число x является представимым.
Все представимые вещественные числа x, для которых |x|≥2^P, суть целые числа. Все представимые вещественные числа x, для которых |x|≥2^P−1, суть целые или полуцелые числа. Все представимые вещественные числа x, для которых |x|≥2^P+1, суть четные целые числа. (Если q_max>P+1, то это утверждение можно очевидным образом продолжить: например, все представимые вещественные числа x, для которых |x|≥2^P+2, суть целые числа, кратные 4, и так далее.)

(Возможность представлять полуцелые числа следует из неравенства q_min<0, оговоренного в требовании 1: здесь не обязательно привлекать требование 2.)

Лемма о монотонности сложения

Важным следствием требования 4 к реализации компьютерной арифметики является утверждение, что функция f(x)=a⊕x является монотонной. Или, другими словами, справедлива следующая

Лемма о (нестрогой) монотонности сложения: для любых представимых чисел a, b, c из b≤c следует a⊕b≤a⊕c.

Ради точности следует добавить, что это действительно следует из наших требований только при отсутствии переполнений, т.е. когда обе суммы a+b и a+c не превосходят по модулю максимального представимого вещественного числа. Практические компьтерные реализации вещественных чисел обычно обеспечивают правильность леммы и в случае переполнений, при этом множество представимых чисел дополняется двумя «бесконечными» значениями +∞/−∞ с естественным обобщением на эти «числа» отношения порядка (в частности, это так в языке Java). Однако, мы не будем пользоваться леммой для таких «крайних» ситуаций — во всех случаях, где мы на нее ссылаемся, значения используемых сумм не превосходят или (в худшем случае) близки по модулю к 2^P, а согласно даже ослабленному требованию 3 (q_max≥P) это, как минимум, почти в два раза меньше максимально представимого вещественного числа 2^q_max(1+(2^P−1)/2^P).

Лемма легко следует из того, что a+b≤a+c, а числа a⊕b и a⊕с являются ближайшими представимыми числами к a+b и a+c. Действительно, допустим, a⊕b>a⊕c, и пусть m=((a⊕b) + (a⊕c))/2 есть середина промежутка между этими двумя числами. Всякое число, для которого ближайшим представимым числом является a⊕b, лежит ближе к a⊕b, чем к a⊕c, т.е. справа от m (в нестрогом смысле); поэтому, в частности, a+b≥m. Всякое число, для которого ближайшим представимым числом является a⊕c, лежит ближе к a⊕c, чем к a⊕b, т.е. слева от m (в нестрогом смысле); поэтому, в частности, a+с≤m. Но раз a+b≤a+c, то это возможно лишь в единственном случае a+b=a+c=m, а тогда получилось бы a⊕b=a⊕c. Противоречие. Лемма доказана.

Заметим, что компьютерное сложение является лишь нестрого монотонным. Иначе говоря, вполне может быть, что b<c, однако, в силу ошибок округления, a⊕b=a⊕c.

Теорема округления

Введем следующие обозначения.

Для любого вещественного числа x обозначим ⌊x⌋ целую часть x, т.е. наибольшее целое число, не превосходящее x. Это традиционное математическое обозначение, касающееся математических вещественных чисел.

Для любого представимого вещественного числа x обозначим round(x) величину ⌊x⊕0.5⌋ — результат выполнения операторов, традиционно используемых в языках программирования для выполнения округления к ближайшему целому число (результат, впрочем, вещественный, еще не приведенный к целому типу данных).

Наконец, для любого представимого вещественного числа x обозначим prev(x) предыдущее представимое вещественное число, т.е. наибольшее представимое вещественное число, строго меньшее x. (Для наименьшего из представимых вещественных чисел функция prev(x) не определена.)

Теорема округления. Пусть x — произвольное представимое вещественное число. Если выполнены два условия:

|x|<2^P либо |x|>2^P+1,

x≠prev(0.5),

то справедливо математически точное равенство round(x)=⌊x+0.5⌋; другими словами, ⌊x⊕0.5⌋=⌊x+0.5⌋. В особом случае x=prev(0.5) результатом round(x) может быть либо 0, либо 1.

Для стандартного IEEE-типа double «магическое» число prev(0.5) приблизительно равно 0.4999999999999999444888. Для стандартного IEEE-типа float оно приблизительно равно 0.4999999702.

Доказательство.

Пусть вначале |x|>2^P+1. По лемме о представимости целых и полуцелых, такое x обязательно целое и четное, причем следующее представимое вещественное число не может быть меньше x+2 (мы пользуемся тем, что |x|≠2^P+1). Поэтому единственным ближайшим представимым числом к x+0.5 является само x, т.е. x⊕0.5=x и ⌊x⊕0.5⌋=⌊x⌋=x. Равенство ⌊x+0.5⌋=x очевидно непосредственно.

С этого момента можно считать, что |x|<2^P.

Сразу исключим ситуацию, когда x является целым. Действительно, согласно лемме о представимости целых и полуцелых, все целые и полуцелые числа в диапазоне |y|≤2^P являются представимыми. Поэтому, если x является целым, то точная сумма x+0.5 является представимой (ибо |x+0.5|≤2^P), то есть x⊕0.5=x+0.5: в этом случае доказывать нечего.

Итак, можно считать, что x не является целым.

Обозначим k целую часть x: k=⌊x⌋. Очевидно, −2^P≤k≤2^P−1. Усилим это неравенство до −2^P≤k≤2^P−2: действительно, если k=2^P−1, то по лемме о представимости целых и полуцелых единственное нецелое x, имеющее такую целую часть, есть 2^P−0.5, а тогда x⊕0.5=x+0.5=2^P и доказывать нечего.

Теперь рассмотрим два возможных варианта: x≥k+0.5 или x<k+0.5.

В первом случае (x≥k+0.5) заметим: раз k≤2^P−2, то по лемме о представимости целых и полуцелых все три числа k+0.5, k+1 и k+1.5 суть представимые вещественные числа. Лемма о монотонности сложения позволяет заключить: поскольку k+0.5≤x≤k+1, то k+1≤x⊕0.5≤k+1.5. Это значит, что ⌊x⊕0.5⌋=k+1. Этому же значению равно математическое выражение ⌊x+0.5⌋. Таким образом, в рассматриваемом случае теорема доказана.

Второй случай (x<k+0.5) требует более подробного анализа. Заметим, что теперь мы можем утверждать: оба числа x и x+0.5 суть нецелые числа, принадлежащие одному и тому же единичному интервалу k..k+1: k<x<k+0.5<x+0.5<k+1. Иначе говоря, не существует ни одного целого числа N, лежащего между x и x+0.5, т.е. такого, что x≤N≤x+0.5.

Разберем три возможных ситуации.

I. Пусть |x|≥1. Согласно общей формуле представления вещественных чисел (см. 1-е требование к представлению вещественных чисел в начале статьи) запишем x = ±2^q(1+m/2^P), причем 0≤q<P (поскольку 1≤|x|<2^P). Соответственно,

x+0.5 = ±2^q(1+m/2^P±1/2^q+1) = ±2^q(1+(m±2^P−1−q)/2^P) = ±2^q(1+m'/2^P), где m'=m±2^P−1−q

(знак ± либо везде означает +, либо везде означает −). Если предположить, что m'≤0, то получится, что целое число N=±2^q лежит между нецелыми числами x=±2^q(1+m/2^P) и x+0.5=±2^q(1+m'/2^P), а это противоречит сказанному выше. Если предположить, что m'≥2^P, то получится, что целое число N=±2^q+1 лежит между этими числами, и это тоже противоречит сказанному выше. Остается допустить, что 0<m'<2^P, а это значит, согласно 1-му требованию к представлению вещественных чисел, что число x+0.5 тоже является представимым. Следовательно, x⊕0.5=x+0.5, и в рассматриваемом случае теорема доказана.

Проведенный анализ можно неформально описать следующим образом: поскольку порядок слагаемого x (т.е. q) больше порядка слагаемого 0.5 (т.е. −1), то алгоритм вещественного сложения вначале уменьшит порядок числа 0.5, разделив его мантиссу на соответствующую степень двойки 2^q+1 (попросту сдвинув на q+1 бит), а затем сложит либо вычтет мантиссы. Условие теоремы |x|<2^P позволяет быть уверенным, что при сдвиге мантиссы не произойдет потери точности — в следующем разделе «Практическое выполнение округления» мы увидим, что при |x|>2^P потеря точности может привести к печальным последствиям, т.е. к неверному округлению. Потеря точности в принципе могла бы произойти далее, если при сложении мантисс происходит переполнение или при вычитании мантисс получается отрицательное число, однако такое возможно лишь в ситуации, когда x лежит в правой половине целого единичного интервала, т.е. x≥k+0.5 (k=⌊x⌋). Эту ситуацию мы разобрали выше и доказали, что в этом случае результат округления всегда правильный, даже несмотря на возможную потерю точности.

II. Пусть −1<x≤0, точнее −1<x<−0.5 (вариант x≥⌊x⌋+0.5 мы уже разобрали выше). Согласно общей формуле представления вещественных чисел запишем x = −2⁻¹(1+m/2^P). Соответственно,

x+0.5 = −2⁻¹*m/2^P = −m/2^P+1

Такое число заведомо является представимым — см. 2-е требование к представлению вещественных чисел в начале статьи (q_min<−P) и следствие из него. Таким образом, снова x⊕0.5=x+0.5, и в рассматриваемом случае теорема доказана.

III. Осталось проверить диапазон 0<x<1, точнее 0<x<0.5 (вариант x≥⌊x⌋+0.5 мы уже разобрали выше). Случай x=prev(0.5) здесь оказывается особым. Действительно, согласно общей формуле представления вещественных чисел,

prev(0.5) = 2⁻²(1+(2^P−1)/2^P) = 0.5−1/2^P+2.

Поэтому

prev(0.5)+0.5 = 1−1/2^P+2 = 2⁻¹(1+(2^P−1)/2^P+Δ), где Δ=1/2^P+1,

а такое число уже не является представимым. Оно имеет два ближайших представимых числа: 2⁻¹(1+(2^P−1)/2^P) и 1. Если оно будет округлено к меньшему из двух представимых чисел (т.е. к первому), то получится ⌊prev(0.5)⊕0.5⌋=0. Если оно будет округлено к 1 (в языке программирования Java гарантирован именно такой результат, в большинстве других языков это также наиболее вероятно), то получится ⌊prev(0.5)⊕0.5⌋=1. Очевидно, при этом ⌊prev(0.5)+0.5⌋=0.

Пусть теперь, как и оговорено в условии теоремы, x≠prev(0.5). Иначе говоря, 0<x≤A=prev(prev(0.5)). Легко убедиться, что для числа A=prev(prev(0.5)) имеем

A = 2⁻²(1+(2^P−2)/2^P) = 0.5−1/2^P+1,
A+0.5 = 1−1/2^P+1 = 2⁻¹(1+(2^P−1)/2^P),

т.е. число A+0.5 является представимым и при этом меньшим 1. Значит, A⊕0.5=A+0.5<1. Лемма о монотонности сложения позволяет заключить: поскольку 0≤x≤A, то 0.5≤x⊕0.5≤A⊕0.5<1. Следовательно, ⌊x⊕0.5⌋=0=⌊x+0.5⌋.

Теорема доказана.

Полезно заметить, что на большинстве вычислительных систем, в частности, в языке программирования Java справедлива несколько более сильная форма теоремы округления. А именно, неравенства из первого условия теоремы можно сделать нестрогими:

Усиленная теорема округления. Пусть x — произвольное представимое вещественное число. Если выполнены три условия:

|x|≤2^P либо |x|≥2^P+1,

x≠prev(0.5),

либо x=−2^P, либо x=2^P+1, либо компьютерное сложение реализовано таким образом, что 2^P⊕0.5=2^P и −2^P+1⊕0.5=−2^P+1 (обе суммы являются непредставимыми, согласно лемме о представимости целых и полуцелых),

то справедливо математически точное равенство round(x)=⌊x+0.5⌋; другими словами, ⌊x⊕0.5⌋=⌊x+0.5⌋. В особом случае x=prev(0.5) результатом round(x) может быть либо 0, либо 1.

Свойство сложения, оговоренное в новом (3-м) условии, например, гарантировано в случае, если при округлении полуцелой непредставимой суммы к ближайшему целому числу предпочтение всегда отдается четному числу; такая гарантия есть в языке программирования Java. Если быть максимально точными, то требование 2^P⊕0.5=2^P позволяет сделать нестрогим первое неравенство 1-го условия (|x|≤2^P), а требование −2^P+1⊕0.5=−2^P+1 позволяет сделать нестрогим второе неравенство (|x|≥2^P+1).

Усиленная форма теоремы непосредственно вытекает из уже доказанной теоремы округления. Для вещественного числа x добавились всего 4 новых значения, не предусмотренных исходной теоремой округления: это x=±2^P и x=±2^P+1. В случае x=−2^P лемма о представимости целых и полуцелых гарантирует, что x+0.5 является представимым, следовательно, x⊕0.5=x+0.5. В случае x=2^P+1 рассуждения абсолютно аналогичны доказательству теоремы округления: по лемме о представимости целых и полуцелых следующим представимым числом является 2^P+1+2, поэтому единственное ближайшое представимое число для x+0.5 есть 2^P+1, и x⊕0.5 может быть округлено только к нему. Наконец, для оставшихся двух значений 3-е условие непосредственно постулирует истинность теоремы.

Практическое выполнение округления по алгоритму «прибавить 0.5 и взять целую часть»

С учетом изложенного, рассмотрим, как реально будет вычисляться округленное значение в компьютерной программе, если используется метод «прибавить 0.5 и взять целую часть». А именно, рассмотрим фрагмент кода на языке Java, являющийся реализацией библиотечной функции StrictMath.round(double x) (во всех версиях языка, предшествовавших Java 7):

(long) Math.floor(x + 0.5)

Это — фрагмент программного кода, поэтому здесь знак + означает компьютерное сложение ⊕, а не математическую сумму. Функция Math.floor имеет результат типа double, равный математически точной целой части своего аргумента (максимальное представимое вещественное число, которое не превосходит аргумента функции и при этом является целым). Наконец, оператор (long)a приводит значение a типа double к целочисленному типу long, который в языке Java является 64-битовым целым числом со знаком в диапазоне −2⁶³..2⁶³−1. Итого три действия: сложение, Math.floor, приведение к целому типу.

Аналогичный код, разумеется, можно написать в большинстве языков программирования, включая C, C++, Pascal.

Какие «подводные камни» следует учесть, чтобы правильно понять результат такого вычисления?

Прежде всего заметим, что операция Math.floor (или ее эквивалент из библиотек других языков программирования) не порождает никаких проблем. Для всякого представимого вещественного числа a его целая часть ⌊a⌋ также представима (это очевидным образом следует из наших требований к формату вещественных чисел) и может быть аккуратно вычислена. Все классические языки программирования имеют библиотечную функцию, реализующую эту операцию математически точно, так же, как это делает библиотечный метод Java Math.floor.

Таким образом, результат выполнения операторов “Math.floor(x + 0.5)” (до приведения типа) описывается формулой ⌊x⊕0.5⌋.

Первый «подводный камень» очевиден из теоремы округления: в единственном уникальном случае, когда x=prev(0.5) (максимальное представимое вещественное число, меньшее 0.5, в случае типа double приблизительно 0.4999999999999999444888), на выходе функции Math.floor мы можем получить 1.0, что математически некорректно. В случае языка Java такой результат гарантирован, в других языках очень вероятен, если только не используется нестандартная модель округления.

Второй «подводный камень» связан с ограниченной точностью представления больших чисел double. Он проявляется в случаях, когда разрядность целого типа, к которому должен быть приведено округленное значение, превышает P битов, в частности, в приведенном фрагменте кода на языке Java (Java-тип long 64-разрядный, что означает 63 бита для представления беззнаковых значений, между тем P=52). Если значение x по модулю лежит в диапазоне 2^P≤|x|≤2^P+1, то вполне вероятен (а на Java даже гарантирован) удивительный результат: нечетные целые числа будут «округлены» к следующим за ними четным числам! Действительно, согласно лемме о представимости целых и полуцелых, в этом случае x есть целое число, однако полуцелое число x+0.5 уже непредставимо. Поэтому при вычислении x⊕0.5 произойдет потеря точности, и сумма будет округлена к одному из двух ближайших целых. Скорее всего, им окажется четное число (а вовсе не исходное целое x, если x было нечетным). Это четное целое число будет передано в функцию Math.floor и затем аккуратно возвращено без изменений в качестве результата округления.

Конкретная реализация сложения, отличающаяся от спецификации языка Java, может использовать другие правила округления для случая, когда имеется два представимых числа, одинаково близких к желаемой сумме — например, предпочитать нечетный результат или выбирать направление случайно. Легко видеть, что при большинстве подходов часть значений из диапазона 2^P≤|x|≤2^P+1 будут обработаны неправильно.

Если x отличен от prev(0.5) и не лежит в «опасных» диапазонах 2^P≤x≤2^P+1, −2^P≥x≥−2^P+1, то согласно теореме округления мы имеем ⌊x⊕0.5⌋=⌊x+0.5⌋, и результат операторов “Math.floor(x + 0.5)” соответствует ожиданиям. Следует заметить, что на большинстве систем (включая язык Java) при возможности округлить результат сложения к двум ближайшим представимым числам предпочтение отдается четному числу, и поэтому граничные значения x=±2^P+1 и x=±2^P также вполне безопасны, т.е. дают правильный результат. Более того, усиленная теорема округления гарантирует, что значения x=−2^P и x=2^P+1 безопасны в любом случае.

Третий и последний «подводный камень» связан с приведением вещественного результата к целому типу long. А именно, если число x превосходит по модулю максимально возможное целое число данного типа, то результат может оказаться неожиданным. В языке Java такое приведение сгенерирует ближайшее целое число, представимое типом long, т.е. 2⁶³−1 в случае очень большого положительного x либо −2⁶³ в случае очень большого по абсолютной величине отрицательного x. В других языках поведение может быть иным.

Можно подытожить указанные соображения в виде следующей рекомендации.

Округление вещественного числа к значению целого типа, основанное на Java-коде
(long)Math.floor(x + 0.5)
или эквивалентном коде иного языка, рекомендуется применять только в случае, когда гарантированы три условия:

|x|<2^P или |x|>2^P+1 (для стандартного IEEE-типа double P=52); на практике обычно удобнее всегда требовать выполнения первого неравенства;

результат округления принадлежит к диапазону представимых чисел целого типа используемого языка (в случае языка Java и типа long это является следствием первого неравенства из предыдущего условия);

при этом x≠prev(0.5) (для стандартного IEEE-типа double это приблизительно 0.4999999999999999444888).

В большинстве реализаций, включая язык Java, неравенства первого условия можно ослабить до нестрогих неравенств |x|≤2^P либо |x|≥2^P+1.

Случай x=prev(0.5), разумеется, при необходимости всегда можно проверить отдельно и округлить данное уникальное значение к 0, а не к 1. На языке Java, начиная с версии Java 6, для получения числа prev(0.5) можно воспользоваться библиотечным вызовом “Math.nextAfter(0.5, 0.0)”.

Теорема вычитания

Приведем в качестве примера одно полезное следствие теоремы округления, которое реально использовалось автором в некоторых прикладных алгоритмах.

Теорема вычитания. Пусть N=2^P, a и b — произвольные представимые вещественные число, такие, что |a|<N, |b|<N и |a−b|≤N. Тогда |round(a)−round(b)|≤N; другими словами, |⌊a⊕0.5⌋−⌊b⊕0.5⌋|≤2^P.

(Здесь подразумевается математически точная разность, а не результат компьютерного вычитания. Как и ранее, P есть число значимых бит в мантиссе — см. 1-е требование к представлению вещественных чисел в начале статьи.)

Прежде всего заметим, что у этой теоремы есть довольно очевидный математический аналог, оперирующий чисто математическими операциями. А именно:

«Математическая» теорема вычитания. Пусть a и b — произвольные вещественные число, N — произвольное неотрицательное целое число и |a−b|≤N. Тогда |⌊a+0.5⌋−⌊b+0.5⌋|≤N.

Докажем это. Пусть для определенности a≥b. Положим a'=b+N; по условию теоремы a≤a'. У чисел a' и b одинаковая дробная часть Δ, а целые части отличаются ровно на N. Если Δ<0.5, то прибавление 0.5 не изменит целую часть чисел a' и b, а если Δ≥0.5, то прибавление 0.5 увеличит целую часть обоих чисел на единицу. В обоих случах получается, что ⌊a'+0.5⌋−⌊b+0.5⌋=N. Остается заметить, что функция f(x)=⌊x⌋ является (нестрого) монотонной, поэтому из a≥b и ⌊a'+0.5⌋=⌊b+0.5⌋+N вытекает ⌊b+0.5⌋≤⌊a+0.5⌋≤⌊b+0.5⌋+N, что и доказывает теорему.

От одной этой теоремы, однако, было бы довольно мало проку при практических вычислениях, где результаты могут отличаться от «идеальных» математических. К счастью, при N=2^P (на самом деле и при некоторых меньших значениях N) справедлива сформулированная перед этим «компьютерная» формулировка той же теоремы. Перейдем к доказательству теоремы вычитания в исходной, «компьютерной» формулировке.

Пусть вновь для определенности a≥b. В силу леммы о монотонности сложения имеем также a⊕0.5≥b⊕0.5, а в силу монотонности функции f(x)=⌊x⌋ также ⌊a⊕0.5⌋≥⌊b⊕0.5⌋, так что нам нужно доказать, что ⌊a⊕0.5⌋−⌊b⊕0.5⌋≤2^P.

Заметим, что в силу леммы о представимости целых и полуцелых неравенства |a|<2^P и |b|<2^P можно переписать так: |a|≤2^P−0.5 и |b|<2^P−0.5 (таково наибольшее представимое число, меньшее 2^P). По лемме о монотонности сложения, отсюда следует, что a⊕0.5≤2^P, b⊕0.5≥−2^P+1, а это значит, что верны также неравенства ⌊a⊕0.5⌋≤2^P, ⌊b⊕0.5⌋≥−2^P+1.

Если ⌊a⊕0.5⌋=⌊a+0.5⌋ и ⌊b⊕0.5⌋=⌊b+0.5⌋, то достаточно сослаться на только что доказанную «математическую» теорему вычитания. А в силу ограничений |a|<2^P и |b|<2^P теорема округления позволяет утверждать, что эти равенства могут нарушиться только в случаях, когда a=prev(0.5) или b=prev(0.5).

Если b=prev(0.5), то по теореме округления ⌊b⊕0.5⌋ может быть равно либо 0, либо 1, т.е. в любом случае неотрицательно. Отсюда сразу следует, что ⌊a⊕0.5⌋−⌊b⊕0.5⌋≤⌊a⊕0.5⌋≤2^P (ибо, как показано выше, ⌊a⊕0.5⌋≤2^P).

Если a=prev(0.5), то по теореме округления ⌊a⊕0.5⌋ может быть равно либо 0, либо 1, т.е. в любом случае не превосходит 1. Отсюда сразу следует, что ⌊a⊕0.5⌋−⌊b⊕0.5⌋≤1−⌊b⊕0.5⌋≤2^P (ибо, как показано выше, ⌊b⊕0.5⌋≥−2^P+1).

Теорема доказана.

Полезно заметить, что на большинстве вычислительных систем, в частности, в языке программирования Java справедлива несколько более сильная форма «компьютерной» теоремы вычитания. А именно, можно разрешить числам a и b быть равными крайним значениям диапазна ±2^P:

Усиленная теорема вычитания. Пусть N=2^P, a и b — произвольные представимые вещественные число, такие, что |a|≤N, |b|≤N и |a−b|≤N. Пусть при этом компьютерное сложение реализовано таким образом, что 2^P⊕0.5=2^P (эта сумма являются непредставимой, согласно лемме о представимости целых и полуцелых). Тогда |round(a)−round(b)|≤N; другими словами, |⌊a⊕0.5⌋−⌊b⊕0.5⌋|≤2^P.

Требуемое свойство компьютерного сложения, например, гарантировано в случае, если при округлении полуцелой непредставимой суммы к ближайшему целому числу предпочтение всегда отдается четному числу; такая гарантия есть в языке программирования Java.

Докажем эту теорему. Как и раньше, положим a≥b. Как и раньше, мы можем утверждать, что ⌊a⊕0.5⌋≤2^P, поскольку для единственного возможного нового значения a=2^P, добавившегося в новой формулироыке теоремы, сама формулировка оговаривает, что a⊕0.5=2^P. Однако, для b лемма о монотонности сложения уже автоматически не даст «удобного» неравенства ⌊b⊕0.5⌋≥−2^P+1, которое мы имели раньше.

Как и раньше, нам достаточно разобрать случаи, когда не выполнено хотя бы одно из равенств ⌊a⊕0.5⌋=⌊a+0.5⌋ и ⌊b⊕0.5⌋=⌊b+0.5⌋ (иначе можно сослаться на «математическую» теорему вычитания). Благодаря усиленной теореме округления мы вновь можем утверждать, что это может случиться только в случаях, когда a=prev(0.5) или b=prev(0.5).

Случай b=prev(0.5) разбирается в точности так же, как и раньше, поскольку у нас сохранилось неравенство ⌊a⊕0.5⌋≤2^P.

В случае a=prev(0.5) заметим, что в силу исходного точного неравенства |a−b|≤2^P мы имеем b≥−2^P+prev(0.5). По лемме о представимости целых и полуцелых ближайшее представимое число, большее −2^P, есть −2^P+0.5, так что можно утверждать b≥−2^P+0.5. А отсюда по лемме о монотонности сложения получается b⊕0.5≥−2^P+1, следовательно, ⌊b⊕0.5⌋≥−2^P+1. Далее рассуждаем, как в прошлом доказательстве. По усиленной теореме округления ⌊a⊕0.5⌋ может быть равно либо 0, либо 1, т.е. в любом случае не превосходит 1. Следовательно, ⌊a⊕0.5⌋−⌊b⊕0.5⌋≤1−⌊b⊕0.5⌋≤2^P.

Теорема доказана.

Учет теоремы округления в Java 7

Разработчики языка программирования Java учли теорему округления в очередной версии языка Java 7. В предыдущих версиях вплоть до Java 6 стандартные функции Math.round и StrictMath.round были реализованы так, как описано в данной статье, т.е. с помощью операторов

(long) Math.floor(x + 0.5)

для типа double или

(int) Math.floor(x + 0.5f)

для типа float. Однако в версии языка Java 7 эти стандартные функции были изменены, равно как и документация к ним (что вряд ли можно считать идеальным решением с точки зрения совместимости). А именно, в Java 7 эти функции для типов данных double и float реализованы следующим образом:

    public static long round(double a) {
        if (a != 0x1.fffffffffffffp-2) // greatest double value less than 0.5
            return (long)floor(a + 0.5d);
        else
            return 0;
    }

    public static int round(float a) {
        if (a != 0x1.fffffep-2f) // greatest float value less than 0.5
            return (int)floor(a + 0.5f);
        else
            return 0;
    }

Иначе говоря, особый случай prev(0.5) в языке Java отныне обрабатывается корректно: возвращается математически точное значение ⌊x+0.5⌋, равное 0. Впрочем, эти функции по-прежнему работают неверно (не в соответствии с естественными ожиданиями) в «опасных» диапазонах 2^P≤x≤2^P+1, −2^P≥x≥−2^P+1, упомянутых в теореме округления: здесь нечетные целые числа «округляются» к соседним четным. Быть может, в следующих версиях языка Java это поведение тоже будет скорректировано.

Ноябрь 2012 г.

AlgART / Тексты