Раздел 23 · Rust

Вещественные числа и IEEE-754

middle-senior~30 мин

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

Вещественные числа и IEEE-754

0.1 плюс 0.2 не равно 0.3, и это не баг, а следствие формата. Знак, экспонента, мантисса, денормали, NaN и бесконечности: разбираем IEEE-754 по битам и выводим правила, когда float уместен, а когда нужен целочисленный расчёт.

Идея

Классика, с которой начинается любой разговор про вещественные числа:

let sum = 0.1 + 0.2;
println!("{sum}");            // 0.30000000000000004
assert!(sum != 0.3);          // проходит

В JavaScript ты списывал это на «такой уж язык». Rust лишает отговорки: здесь строгие типы, никакой автоконвертации, а результат тот же. Потому что дело не в языке, а в формате IEEE-754, общем для всех языков и процессоров.

Сам формат описывается одной фразой: это научная запись в двоичной системе с фиксированным числом разрядов. Та самая запись из физики, 6.022 умножить на 10 в 23-й, только основание не 10, а 2, и под мантиссу с экспонентой отведено жёсткое число бит. Из этой фразы выводится всё поведение float, включая знаменитую сумму выше: 0.1 в двоичной системе это бесконечная периодическая дробь, как 1/3 в десятичной, и хранить её точно невозможно в принципе, сколько бит ни выдай.

Формат не врёт, он округляет. Урок про то, где именно: после него ты сможешь сказать, какие числа float хранит точно, какие нет, почему сравнение на равенство опасно и что делать с сортировкой, в которую затесался NaN. После целых из прошлого урока это вторая, последняя половина числового мира.

Знак, экспонента, мантисса

Раскладка f32 по битам: один бит знака, восемь бит экспоненты, двадцать три бита мантиссы.

s eeeeeeee mmmmmmmmmmmmmmmmmmmmmmm
1    8                23           = 32 бита

Значение собирается по формуле: минус единица в степени s, умножить на 1.m, умножить на 2 в степени e минус 127. Три детали формулы стоят разворота.

Экспонента хранится со смещением: в поле лежит не сама степень, а степень плюс 127. Так сделано, чтобы порядок float как битовых строк почти совпадал с числовым порядком; реальные степени для f32 бегут от минус 126 до плюс 127.

Мантисса хранится без целой части: в двоичной научной записи нормализованное число всегда начинается с единицы, 1.что-то, поэтому единицу не хранят, она подразумевается. Это скрытый двадцать четвёртый бит точности, бесплатный.

И знак это отдельный бит, а не дополнительный код: float кодирует знак как у людей, знак и величина. Отсюда, забегая вперёд, существование минус нуля.

В Rust на всё это можно посмотреть без единого unsafe, инструментами позапрошлого урока:

let x = 1.0f32;
let bits = x.to_bits(); // тот же u32, та же память, другая интерпретация

println!("{bits:#034b}"); // 0b00111111100000000000000000000000

let sign = bits >> 31;
let exponent = (bits >> 23) & 0xFF;
let mantissa = bits & 0x7F_FFFF;

assert_eq!(sign, 0);
assert_eq!(exponent, 127); // 127 - 127 = 0, степень нулевая
assert_eq!(mantissa, 0);   // 1.0 без дробной части

Единица это плюс, степень ноль, мантисса 1.0: всё сходится. to_bits и from_bits будут нашим микроскопом весь урок, а маски и сдвиги ты тренировал два урока назад. Теперь посмотрим на знаменитую десятую долю:

let x = 0.1f32;
println!("{:#x}", x.to_bits());  // 0x3dcccccd
println!("{x:.20}");             // 0.10000000149011611938

Мантисса ccccc... с хвостом d: период 1100 обрезан на двадцать третьем бите и округлён вверх. В переменной лежит не 0.1, а ближайшее представимое число, на полтора десятимиллионных процента больше. Дальше арифметика честно работает с этим соседом, и ошибки двух таких соседей в сумме 0.1 и 0.2 складываются ровно в тот хвостик, который ты видел в начале урока.

f64 устроен так же, только бит больше: знак, одиннадцать бит экспоненты со смещением 1023, пятьдесят два бита мантиссы. Точность около шестнадцати десятичных цифр против семи у f32. Выбор между ними прост: по умолчанию f64, как и делает Rust для литерала 3.14; f32 это осознанная экономия памяти и пропускной способности в графике, машинном обучении и больших массивах.

Специальные значения

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

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

let inf = f32::MAX * 2.0;
assert_eq!(inf, f32::INFINITY);
assert_eq!(1.0f32 / 0.0, f32::INFINITY);
assert_eq!(-1.0f32 / 0.0, f32::NEG_INFINITY);

Заметь контраст с целыми: float при переполнении не паникует и не заворачивается, а уезжает в бесконечность. Семантика «прилип к краю» из прошлого урока, только край здесь не максимум типа, а отдельное значение.

Экспонента из всех единиц с ненулевой мантиссой это NaN. Он рождается там, где у операции нет осмысленного ответа, и обладает двумя свойствами, из-за которых о нём нельзя не знать. Первое: NaN заразен, любая арифметика с ним даёт NaN, и ошибка тихо проедет через весь расчёт до самого вывода. Второе: NaN не равен ничему, включая самого себя.

let nan = f32::NAN;
assert!(nan != nan);          // единственное значение в языке с таким свойством
assert!(!(nan < 1.0) && !(nan > 1.0) && nan != 1.0); // несравним ни с чем
assert!(nan.is_nan());        // проверка только так

Сравнение с NaN всегда ложь, поэтому проверять на него можно только методом is_nan. Это не каприз Rust, а требование стандарта IEEE-754, одинаковое во всех языках, и через раздел оно объяснит нам, почему float так неудобно сортировать.

Экспонента из всех нулей это два жителя. При нулевой мантиссе ноль, причём из-за знакового бита их два: 0.0 и -0.0. Они равны при сравнении, но различимы по битам и по поведению, единица на минус ноль даёт минус бесконечность. При ненулевой мантиссе денормали: совсем крошечные числа, у которых скрытая единица выключена. Без них между нулём и самым маленьким нормализованным числом зияла бы дыра, и вычитание двух разных маленьких чисел могло бы дать ровно ноль; денормали заполняют дыру, жертвуя точностью постепенно. Для прикладного кода важно одно: появление денормалей в расчёте означает, что ты работаешь на самом дне диапазона и точность уже расползлась.

Точность и сравнение

Вернёмся из зоопарка к главному практическому вопросу: насколько float точен и как с этим жить.

Ключевая картинка такая: представимые числа разбросаны по оси неравномерно. Мантисса даёт фиксированные 24 бита значащих цифр для f32, а экспонента двигает их по оси, поэтому между 1 и 2 представимых чисел столько же, сколько между 1024 и 2048: шаг сетки растёт вместе с числом. Около единицы шаг f32 примерно одна десятимиллионная, а после 16 777 216, это 2 в степени 24, шаг становится больше единицы, и целые числа начинают проскакивать:

let big = 16_777_216f32;          // 2^24
assert_eq!(big + 1.0, big);       // единица меньше шага сетки, прибавка исчезла

У f64 тот же порог это 2 в степени 53, около девяти квадриллионов. Это число ты уже встречал: Number.MAX_SAFE_INTEGER в JavaScript, где все числа и есть f64. Теперь ты знаешь, откуда оно берётся: это не предел типа, а граница, после которой шаг сетки превышает единицу и целые перестают представляться точно.

Из шага сетки следует и правило сравнения. Два математически равных выражения, вычисленные разными путями, приходят к соседним узлам сетки, поэтому == для результатов расчёта почти всегда ошибка. Сравнивай с допуском:

fn approx_eq(a: f64, b: f64, eps: f64) -> bool {
    (a - b).abs() <= eps * a.abs().max(b.abs()).max(1.0)
}

assert!(approx_eq(0.1 + 0.2, 0.3, 1e-12));

Допуск масштабируется по величине операндов: сетка ведь тоже масштабируется, и абсолютный допуск, годный около единицы, бессмыслен около миллиарда. Константа f64::EPSILON, шаг сетки около единицы, служит ориентиром для выбора eps, но не готовым ответом: после цепочки операций накопленная ошибка обычно больше одного шага.

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

let a = 1.000_000_1f32;
let b = 1.0f32;
println!("{:e}", a - b); // 1.1920929e-7, а математически 1e-7: ошибка 19 процентов

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

Отсюда же главное табу: float не для денег. Сумма 0.1 + 0.2 уже не равна 0.3, а бухгалтерия требует, чтобы копейки сходились до бита, и проверяющему не объяснишь про сетку представимых чисел. Деньги считают в целых минимальных единицах, копейках или центах, инструментами прошлого урока: u64 и checked_-семейство. Когда нужны дробные доли процента и десятичная точность, берут десятичные типы вроде крейта rust_decimal, внутри которых целая мантисса и десятичная экспонента. Тот же принцип распространяется на блокчейн, и там он жёстче: консенсус требует бит-в-бит одинакового результата на всех узлах, поэтому балансы в эфире это целые wei, 10 в восемнадцатой степени на один эфир, а float в коде контрактов не существует как класс.

Где float уместен и незаменим: физика, геометрия, графика, сигналы, машинное обучение, всё, где величины по природе непрерывны и измерены с погрешностью, на фоне которой ошибка округления шум. Диапазон f64 до 10 в 308-й при шестнадцати значащих цифрах целым типам недоступен в принципе.

Тотальный порядок

Последний сюжет связывает float с системой трейтов из блока про язык и объясняет ежедневное неудобство.

Попробуй отсортировать вектор float очевидным способом:

let mut xs = vec![3.0f64, 1.0, 2.0];
xs.sort(); // не компилируется: f64 не реализует Ord

Ошибка не каприз. Ord обещает тотальный порядок: любые два значения сравнимы. NaN это обещание ломает, он несравним ни с чем, поэтому f64 честно реализует только PartialOrd. По той же причине у float нет Eq и Hash: рефлексивность равенства разбита тем самым nan != nan, и ключом HashMap или BTreeMap float быть не может. Компилятор не вредничает, а отказывается подписывать контракт, который тип не выполняет; это та же честность типов, что и везде в Rust, просто здесь она впервые мешает.

Выходов два. Когда NaN в данных гарантированно нет, сравнивай через partial_cmp с явным решением, что делать при сюрпризе, или пользуйся sort_by(|a, b| a.partial_cmp(b).unwrap()), осознанно соглашаясь на панику. Но с 1.62 в std есть выход лучше:

let mut xs = vec![3.0f64, f64::NAN, 1.0, 2.0];
xs.sort_by(f64::total_cmp);
// 1.0, 2.0, 3.0, NaN: отсортировалось, NaN ушёл в конец

total_cmp реализует тотальный порядок из самого стандарта IEEE-754: минус NaN, минус бесконечность, числа по порядку, минус ноль перед плюс нулём, плюс бесконечность, плюс NaN. Внутри он сравнивает битовые представления с поправкой на знак, та самая близость порядка бит к порядку чисел, ради которой экспоненту хранят со смещением. От оператора < он отличается ровно в патологиях: NaN получает место в строю, а нули перестают быть равными. Для сортировки, поиска максимума и построения индексов это то, что нужно; для математических сравнений в формулах по-прежнему обычные операторы.

ДЗ

Дальше

Числовой мир закрыт целиком: биты, целые, вещественные. Следующий шаг вниз последний: машинный код. Посмотрим через godbolt, во что компилятор превращает цикл и match, научимся читать дизассемблер, не пугаясь, и узнаем, что осталось от твоих типов после компиляции. Спойлер: ничего, и это лучшая новость урока.