Sannis 28 дек 2017 в 15:17

Обзор реализаций округления в Go

9 мин

34K

Перевод

Привет, Хабр! Меня зовут Олег, я PHP-и-не-только-разработчик в Badoo. Меня часто удивляет, насколько по-разному в языках программирования подходят к составлению стандартной библиотеки. Go — не исключение: отсутствие функции math.Round() меня удивило. Однако, покопавшись в этих ваших интернетах, я выяснил, в чём причина. Этими знаниями я и хотел бы поделиться в своём вольном переводе.

Округление (round) в Go нелегко сделать правильно. Казалось бы, берём float64, отбрасываем дробную часть и прибавляем единицу к получившемуся значению, если дробная часть была >= 0,5. Есть и другие варианты, доступные при поиске “golang round” в Google, многие из которых также можно найти в закрытом тикете на GitHub, где Round() отказались включать в стандартный пакет math. В этой статье я хотел бы изучить эти и другие реализации и проверить их на корректность. Мы увидим, что почти во всех есть ошибки, препятствующие их использованию.

Выбор способа округления

Существует несколько способов округления в зависимости от способа применения результата: округление к меньшему/ большему, округление к меньшему/ большему по модулю, округление к ближайшему целому, округление к ближайшему чётному и т. д… Округление к ближайшему целому, в свою очередь, можно делать по-разному в зависимости от того, какой результат должен получиться, если дробная часть равна 0,5. Я буду рассматривать округление к ближайшему целому, причём 0,5 будет округляться в большую (по модулю) сторону.

Требования к корректной реализации Round() заключаются в следующем:

правильно округляет до ближайшего целого все конечные числа;
поддерживает специальные значения (NaN, Inf, -0), возвращая их без изменений.

Я буду использовать следующие тестовые примеры для проверки корректности, в каждой паре содержатся исходное значение и предполагаемый результат выполнения функции Round():

tests := [][2]float64{
    {-0.49999999999999994, negZero}, // -0.5+epsilon
    {-0.5, -1},
    {-0.5000000000000001, -1}, // -0.5-epsilon
    {0, 0},
    {0.49999999999999994, 0}, // 0.5-epsilon
    {0.5, 1},
    {0.5000000000000001, 1},                         // 0.5+epsilon
    {1.390671161567e-309, 0},                        // denormal
    {2.2517998136852485e+15, 2.251799813685249e+15}, // 1 bit fraction
    {4.503599627370497e+15, 4.503599627370497e+15},  // large integer
    {math.Inf(-1), math.Inf(-1)},
    {math.Inf(1), math.Inf(1)},
    {math.NaN(), math.NaN()},
    {negZero, negZero},
}

В этом списке есть обычные числа, специальные значения и некоторые граничные случаи, с которыми простым алгоритмам сложно справиться. Обратите внимание, что, поскольку мы используем float, мы не можем использовать число 0,49999999999999999 в качестве ближайшего к 0,5, так как из-за ограниченной точности float это число в точности равно 0,5. Вместо этого я использую 0,49999999999999994.

Реализации, предложенные в закрытом тикете, явно не были проверены на подобных данных, часто не работали даже те из них, которые были предложены известными людьми. Это лишний раз доказывает, насколько сложно написать Round().

int(f + 0.5)

Первая реализация, предложенная rsc, выглядела следующим образом:

return int(f + 0.5)

Она некорректно работает с особыми значениями, отрицательными числами, числами больше math.MaxInt64 и числами, близкими к 0,5:

round(-0.5): got: 0, want -1
round(-0.5000000000000001): got: 0, want -1
round(0.49999999999999994): got: 1, want 0
round(4.503599627370497e+15): got: 4.503599627370498e+15, want 4.503599627370497e+15
round(-Inf): got: -9.223372036854776e+18, want -Inf
round(+Inf): got: -9.223372036854776e+18, want +Inf
round(NaN): got: -9.223372036854776e+18, want NaN

Floor() or Ceil()

Второй предложенный вариант учитывал отрицательные числа:

if f < 0 {
    return math.Ceil(f - 0.5)
}
return math.Floor(f + 0.5)

однако продолжал некорректно работать в некоторых случаях:

round(-0.49999999999999994): got: -1, want -0
round(0.49999999999999994): got: 1, want 0
round(4.503599627370497e+15): got: 4.503599627370498e+15, want 4.503599627370497e+15

Первые два теста не проходят, потому что результат разности n — 0,5 равен в точности -1,0, тогда как мы ожидаем получить что-то точно большее, чем -1,0. Если посмотреть на реализацию round в Postgres, можно понять, как решить эту проблему.

Самое интересное, что эта ошибка не является такой уж редкой. До версии 6 точно такая же присутствовала в Java. Хорошо, что с тех пор реализация улучшилась.

int и Copysign

В третьем предложении от minux была предпринята другая попытка решить проблему отрицательных чисел:

return int(f + math.Copysign(0.5, f))

И этот вариант всё равно ломает тесты:

round(-0.49999999999999994): got: -1, want -0
round(0.49999999999999994): got: 1, want 0
round(4.503599627370497e+15): got: 4.503599627370498e+15, want 4.503599627370497e+15
round(-Inf): got: -9.223372036854776e+18, want -Inf
round(+Inf): got: -9.223372036854776e+18, want +Inf
round(NaN): got: -9.223372036854776e+18, want NaN

Как видно, часть тестов стала проходить, однако другие начали падать. Была предпринята попытка улучшить этот алгоритм:

if math.Abs(f) < 0.5 {
        return 0
}
return int(f + math.Copysign(0.5, f))

Однако и она провалилась:

round(4.503599627370497e+15): got: 4.503599627370498e+15, want 4.503599627370497e+15
round(-Inf): got: -9.223372036854776e+18, want -Inf
round(+Inf): got: -9.223372036854776e+18, want +Inf
round(NaN): got: -9.223372036854776e+18, want NaN

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

Мы рассмотрели уже четыре варианта, и в каждом из них нашлись изъяны. Настало время посмотреть, как Round() реализуют авторы различных пакетов.

Kubernetes

Kubernetes 1.7 содержит такую реализацию:

if a < 0 {
    return int32(a - 0.5)
}
return int32(a + 0.5)

Она ломает следующие тесты:

round(-0.49999999999999994): got: -1, want -0
round(0.49999999999999994): got: 1, want 0
round(4.503599627370497e+15): got: 4.503599627370498e+15, want 4.503599627370497e+15
round(-Inf): got: -9.223372036854776e+18, want -Inf
round(+Inf): got: -9.223372036854776e+18, want +Inf
round(NaN): got: -9.223372036854776e+18, want NaN

Судя по тому, что функция возвращает int32, она не предназначена для работы с большими числами. Однако она некорректно работает и с числами, которые близки к 0,5.

Работающие реализации на Go

Round(), используемая в Postgres

Выше я уже упоминал, что в Postgres содержится код функции Round() на C, который работает для всех тестируемых значений. В CockroachDB мы переписали этот код на Go, без комментариев он выглядит следующим образом:

func round(x float64) float64 {
    if math.IsNaN(x) {
        return x
    }
    if x == 0.0 {
        return x
    }
    roundFn := math.Ceil
    if math.Signbit(x) {
        roundFn = math.Floor
    }
    xOrig := x
    x -= math.Copysign(0.5, x)
    if x == 0 || math.Signbit(x) != math.Signbit(xOrig) {
        return math.Copysign(0.0, xOrig)
    }
    if x == xOrig-math.Copysign(1.0, x) {
        return xOrig
    }
    r := roundFn(x)
    if r != x {
        return r
    }
    return roundFn(x*0.5) * 2.0
}

Давайте разберёмся, как он работает. Первые шесть строк обрабатывают особые случаи. Далее мы выбираем roundFn из Ceil и Floor в зависимости от того, положительное число или отрицательное. Далее начинается самое интересное:

x -= math.Copysign(0.5, x)

Этим кодом мы сдвигаем x ближе к нулю.

if x == 0 || math.Signbit(x) != math.Signbit(xOrig) {
    return math.Copysign(0.0, xOrig)
}

Далее мы проверяем, не стал ли x в точности нулём и не поменялся ли у него знак. Это означает что исходное число <= 0,5, в этом случае мы возвращаем ноль с нужным знаком.

if x == xOrig-math.Copysign(1.0, x) {
    return xOrig
}

Эта проверка нужна для очень больших чисел, для которых x-0,5 == x-1,0, в этих случаях мы можем вернуть число неизменённым.

r := roundFn(x)
if r != x {
    return r
}

Далее мы округляем число с помощью Floor() или Ceil() и возвращаем это значение, если оно отличается от x, что может случиться, только если дробная часть входного значения не равна в точности 0,5, так как выше мы вычли 0,5 из него.

return roundFn(x*0.5) * 2.0

Теперь мы знаем, что дробная часть равна 0,5, поэтому нам нужно округлить до ближайшего чётного числа (реализация Round() в Postgres в этом месте отличается от приведённых выше вариантов). Комментарий в коде лучше это описывает:

Dividing input+0.5 by 2, taking the floor and multiplying by 2 yields the closest even number. This part assumes that division by 2 is exact, which should be OK because underflow is impossible here: x is an integer.

Чтобы сохранить оригинальное поведение, этот код можно заменить на следующий:

return xOrig + math.Copysign(0.5, xOrig)

github.com/montanaflynn/stats

Ещё одна работающая реализация содержится в пакете github.com/montanaflynn/stats. Без комментариев она выглядит следующим образом:

func round(input float64) {
    if math.IsNaN(input) {
        return math.NaN()
    }
    sign := 1.0
    if input < 0 {
        sign = -1
        input *= -1
    }
    _, decimal := math.Modf(input)
    var rounded float64
    if decimal >= 0.5 {
        rounded = math.Ceil(input)
    } else {
        rounded = math.Floor(input)
    }
    return rounded * sign
}

Ключевое отличие от предыдущих решений заключается в использовании функции Modf(), которая корректно разделяет целую и дробную части чисел.

Round() в Go 1.10

Спустя несколько месяцев после выхода Go 1.8 появился очередной тикет с просьбой добавить math.Round в Go. В комментариях к нему продолжали появляться некорректно работающие реализации, их число увеличилось до восьми. К счастью, команда Go согласилась добавить math.Round в Go 1.10! И даже появилась работающая реализация:

func Round(x float64) float64 {
    const (
        mask     = 0x7FF
        shift    = 64 - 11 - 1
        bias     = 1023

        signMask = 1 << 63
        fracMask = (1 << shift) - 1
        halfMask = 1 << (shift - 1)
        one      = bias << shift
    )

    bits := math.Float64bits(x)
    e := uint(bits>>shift) & mask
    switch {
    case e < bias:
        // Round abs(x)<1 including denormals.
        bits &= signMask // +-0
        if e == bias-1 {
            bits |= one // +-1
        }
    case e < bias+shift:
        // Round any abs(x)>=1 containing a fractional component [0,1).
        e -= bias
        bits += halfMask >> e
        bits &^= fracMask >> e
    }
    return math.Float64frombits(bits)
}

Для тех, кто не знаком с устройством float (я в их числе), этот код выглядит совершенно непонятно. Попробуем разобраться, что же он делает:

bits := math.Float64bits(x)
e := uint(bits>>shift) & mask

Похоже, что мы берём битовое представление числа, сдвигаем его и применяем маску. Согласно спецификации IEEE 754:

The encoding scheme for these binary interchange formats is the same as that of IEEE 754-1985: a sign bit, followed by w exponent bits that describe the exponent offset by a bias, and p−1 bits that describe the significand.

Рассматривая приведённые выше константы, мы видим, что сдвиг составляет 64 — 11 — 1, что означает 64 бита на число, 11 из которых используются для показателя степени, один — для знака и 52 оставшихся бита — для мантиссы. Это означает, что используемый сдвиг удаляет биты мантиссы, а маска удаляет бит знака, оставляя нас только с показателем степени.

switch {
case e < bias:

В полученном числе показатель степени записан не как он есть, а с прибавлением числа 1023 (это делается для того чтобы записывать отрицательные показатели для очень маленьких чисел), что означает, что мы должны вычесть 1023 из e, вычисленного выше, чтобы получить фактический показатель. Иными словами, если e < bias, то мы имеем отрицательный показатель степени, что означает, что абсолютное значение float должно быть < 1. Действительно, дальше мы видим:

// Round abs(x)<1 including denormals.
bits &= signMask // +-0
if e == bias-1 {
    bits |= one // +-1
}

Здесь бит маскируется знаковым битом, это используется только для сохранения правильного знака: теперь мы можем полностью игнорировать мантиссу. Мы можем это сделать, потому что в этом случае нас интересует только показатель степени. Так как используется основание степени 2, а e < bias, мы знаем, что наименьший показатель, который может быть, равен -1, а 2 ^ -1 = 0,5. Кроме того, мантисса имеет некоторое значение 1.X. Таким образом, в зависимости от показателя наше число находится либо в диапазоне (0,5, 1), либо в диапазоне (0, 0,5). Поэтому во втором случае для правильного округления нам нужно добавить к числу единицу. Фух. Подробнее это описано в википедии.

Теперь разберём второй случай:

case e < bias+shift:

Наверное, вы думаете, что условие в этой ветке должно быть e > bias, чтобы покрыть все случаи с положительным показателем степени. Но вместо этого тут используется только их часть. Использование сдвига здесь особенно интересно, потому что кажется, что оно несравнимо с bias. Первое — это число битов смещения, а второе — численное смещение. Но, поскольку числа с плавающей точкой представлены как (1.мантисса) * 2 ^ X, то если X больше числа битов в мантиссе, мы гарантированно получим значение без дробной части. То есть показатель степени сместил десятичную точку вправо настолько, что мантисса окончательно пропала. Таким образом, выражение в этой ветке игнорирует числа с плавающей точкой, которые уже округлены.

// Round any abs(x)>=1 containing a fractional component [0,1).
e -= bias
bits += halfMask >> e
bits &^= fracMask >> e

Первая строка тут простая: вычитаем bias из e и получаем реальное значение показателя степени. Вторая строка добавляет к значению 0,5. Это работает, потому что старший бит мантиссы добавляет 0,5 к финальной сумме (см. представление в статье “Википедии” ниже). В этом случае эта сумма переполняет 52-битные границы мантиссы, показатель степени будет увеличен на 1. Значение показателя степени не сможет переполниться до знакового бита, так как оно не может быть больше bias+shift из примера выше. В любом случае, дробная часть очищается. Таким образом, если дробная часть была больше или равна 0,5, она будет увеличена на 1, в противном случае будет отброшена. Хитро и не очевидно до тех пор, пока мы не посмотрим глубже.

Заключение

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

Думаю, команда Go приняла правильное решение, добавив функцию Round() в стандартную библиотеку. Без этого мы бы продолжали пользоваться различными некорректными реализациями.

Надеюсь, теперь вам стало ясно, что при работе с float есть много подводных камней, про которые порой забывают даже эксперты. Легко придумать или скопировать откуда-то однострочную реализацию, но сложно написать действительно корректную. Неудивительно, что корректно работающее округление появилось лишь в шестой мажорной версии Java (через 15 лет, прошедших с релиза Java 1.0 до выхода Java 7), и я рад, что Go прошёл этот путь быстрее.

Теги:

Хабы: