Арифметика компьютера: точность, ошибки и цена вычислений¶

Точность, ошибки и инженерные компромиссы в цифровом счёте

Прежде чем начать: три вопроса без ответа¶

Запишите ответы на листке — не обсуждая с соседом.

В базовом курсе мы увидели: 0.1 + 0.2 != 0.3 в Python. Вы знаете симптом. А можете ли вы объяснить, какой именно бит отличает 0.1 + 0.2 от 0.3 в памяти процессора?
int в Python никогда не переполняется. В Java int переполняется по модулю. В C переполнение signed int — неопределённое поведение. Три языка, три ответа. Почему?
Ракета Ariane 5 стоимостью 370 миллионов долларов взорвалась через 37 секунд из-за одной операции с числом. Что именно пошло не так?

На первый вопрос у большинства программистов нет правильного ответа. На второй — удивляет всех, кто привык к одному языку. Третий изменил историю ракетостроения.

От симптома к механизму¶

В базовом курсе мы увидели, что 0.1 + 0.2 ≠ 0.3, что числа в компьютере конечны, и что math.isclose() спасает от наивного сравнения. Это были симптомы.

Эта глава — диагноз и инженерное лечение. Мы разберём, как именно бит кодирует число: от транзистора до регистра процессора. Увидим, где рождается ошибка — не «где-то в компьютере», а в конкретном разряде мантиссы. И научимся принимать инженерное решение: какой тип данных выбрать под конкретную задачу, как оценить погрешность до запуска программы, и когда точность вычисления уже не имеет значения — потому что ограничена точностью измерения.

Предыдущая глава ответила на вопрос «что такое число». Эта глава отвечает на другой: «что компьютер с ним делает — и чего это стоит».

Часть I. Целые числа — физика битов¶

§ 1. Целые числа в памяти: точная арифметика с конечным диапазоном¶

Провокация¶

Напишем на доске: int x = 2147483647; x = x + 1; — каков результат?

Ожидаемый ответ: 2147483648. В Java на типе int получится -2147483648: гарантированное переполнение по модулю $2^{32}$. В C для signed int такое выражение формально приводит к undefined behavior — стандарт не гарантирует результат, хотя на большинстве платформ вы увидите то же значение. Это не баг — это следствие разных контрактов языка с программистом.

§ 1.1 Почему n бит — это ровно 2ⁿ состояний¶

Компьютерная память построена на транзисторах. Каждый транзистор имеет два устойчивых состояния: есть ток / нет тока, заряжен / разряжен. Одно состояние — один бит: 0 или 1. Это физика полупроводников, не соглашение. При n битах возможно ровно $2^n$ различных комбинаций. Никак не больше и не меньше. Это означает: целое число в n-битной ячейке памяти может принять ровно $2^n$ различных значений. Какие именно — определяет выбор кодирования.

§ 1.2 Беззнаковые целые: чистая арифметика по модулю¶

uint8: 8 бит → $2^8$ = 256 значений → диапазон 0..255

uint16: 16 бит → $2^{16}$ = 65 536 значений → диапазон 0..65 535

uint32: 32 бита → $2^{32}$ ≈ 4.3 млрд значений → диапазон 0..4 294 967 295

uint64: 64 бита → $2^{64}$ ≈ 1.8 · $10^{19}$ значений

Все арифметические операции выполняются по модулю $2^n$. 255 + 1 = 0 для uint8 — не ошибка, а математически корректный результат в ℤ/256ℤ.

Арифметика по модулю — это не упрощение и не ошибка языка программирования. Это прямое следствие физики: n транзисторов не могут хранить больше $2^n$ состояний.

Арифметика uint8 в действии¶

200 + 100 = 300 → 300 mod 256 = 44 (переполнение!)

255 + 1 = 256 → 256 mod 256 = 0 (оборот к нулю)

0 - 1 = -1 → -1 mod 256 = 255 (оборот к максимуму)

На Python: import numpy as np; x = np.uint8(200); x + np.uint8(100) → 44

Практическое применение: IP-адрес — это uint32 (четыре байта, каждый 0–255). Именно поэтому IPv4-адресов всего ~4.3 миллиарда — в точности $2^{32}$. Именно поэтому интернет «закончился» в 2011 году и пришлось переходить на IPv6 с uint128.

§ 1.3 Знаковые целые: дополнение до двух¶

Как закодировать отрицательные числа? Первый вариант: один бит — знак, остальные — значение (sign-and-magnitude). Проблема: два нуля (+0 и -0) и разные схемы сложения для положительных и отрицательных. Решение, принятое во всех современных процессорах: дополнение до двух (two's complement). Оно устраняет обе проблемы.

Кодирование в дополнение до двух¶

Кодирование в дополнение до двух. Неотрицательные числа 0..$2^{n-1}-1$: кодируются обычно, в двоичном представлении.

Отрицательное число -x (где x > 0): кодируется как $2^n$ - x.

Практическое правило: инвертировать все биты и прибавить 1.

Диапазон int8: от -128 до +127. Почему асимметрия? -128 = 10000000₂, и у него нет положительного двойника.

Ключевое свойство: сложение в дополнении до двух работает как сложение по модулю $2^n$ — одна схема для всех случаев.

Пример: кодирование −42 в int8¶

Кодирование −42 в int8. 1. Записываем 42 в двоичном: 00101010₂

Инвертируем все биты: 11010101₂
Прибавляем 1: 11010110₂

Проверка: 11010110₂ = 128+64+16+4+2 = 214 в uint8. И действительно: -42 + 256 = 214.

Обратное декодирование: старший бит = 1 → число отрицательное. Значение: -(~11010110₂ + 1) = -00101010₂ = -42.

§ 1.4 Битовые операции как быстрая арифметика¶

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

Битовые операции и их арифметический смысл¶

x << k (сдвиг влево на k): для беззнаковых или неотрицательных целых — умножение на $2^k$.

x >> k (сдвиг вправо): для беззнаковых или неотрицательных — целочисленное деление на $2^k$ с отбросом дробной части. (Для отрицательных знаковых целых поведение зависит от языка.)

x & ($2^k$-1): в типичной двоичной арифметике для беззнаковых или неотрицательных целых даёт остаток от деления на $2^k$. Пример: x & 7 = x mod 8.

x & (x-1): обнуляет младший установленный бит. Применение: подсчёт битов, проверка степеней двойки.

Эти операции — основа хеш-функций, алгоритмов сжатия, криптографических примитивов.

§ 1.5 Стандартные типы и их реальные размеры¶

Тип	Размер	Диапазон	Применение
int8 / uint8	1 байт	−128..127 / 0..255	Байты, пиксели, ASCII
int16 / uint16	2 байта	−32768..32767 / 0..65535	Порты TCP/UDP, аудио
int32 / uint32	4 байта	±2.1·10⁹ / 0..4.3·10⁹	Стандартный int (C/Java)
int64 / uint64	8 байт	±9.2·10¹⁸ / 0..1.8·10¹⁹	long, файловые позиции
Python int	переменный	неограничен	BigInt (подробнее в §2)

Ariane 5: 370 миллионов долларов за 37 секунд¶

4 июня 1996 года. Ariane 5, первый полёт нового европейского носителя. 37 секунд после старта — ракета разрушила себя сама.

Причина: модуль инерциальной навигационной системы (IRS), унаследованный от Ariane 4, содержал код, преобразующий 64-битное число с плавающей точкой (горизонтальная скорость) в 16-битное знаковое целое (int16). Для Ariane 4 это значение никогда не превышало 32 767. Для Ariane 5 с её другой траекторией — превысило.

Преобразование float64 → int16 выбросило исключение переполнения. Исключение не было обработано. Система перешла в диагностический режим и начала передавать диагностические данные вместо полётных. Бортовой компьютер интерпретировал диагностику как данные об ориентации, дал коррекционную команду — ракета резко отклонилась и разрушилась.

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

Типичная ошибка¶

Знаковое/беззнаковое сравнение в C: if (x < y) при x = -1 (int, знаковое) и y = 1 (unsigned int, беззнаковое) даёт FALSE. Компилятор C преобразует -1 к типу unsigned int: -1 в дополнении до двух = 0xFFFFFFFF = 4 294 967 295. Это огромное положительное число, больше 1.

Этот класс ошибок (integer sign confusion) стабильно входит в топ критических уязвимостей безопасности (CWE Top 25).

Код: переполнение и дополнение до двух¶

Что демонстрирует: кодирование отрицательных чисел в two's complement, наблюдение переполнения через numpy.

def twos_complement(n, bits):
    if n < 0:
        n = (1 << bits) + n
    return format(n & ((1 << bits) - 1), f'0{bits}b')
print(twos_complement(-1, 8))    # 11111111
print(twos_complement(-128, 8))  # 10000000
print(twos_complement(127, 8))   # 01111111
# Реальное переполнение через numpy
import numpy as np
x = np.int8(127)
print(x + np.int8(1))            # -128 (переполнение!)
print(np.uint8(255) + np.uint8(1))  # 0 (оборот)

Задачи к §1¶

Базовый. Закодируйте вручную числа -1, -127, -128 в 8-битное дополнение до двух. Проверьте через Python.
Аналитический. Напишите функцию decode_twos(bits_str) → int, которая декодирует строку битов обратно в целое.
Аналитический. Объясните: почему int8(-128) / int8(-1) является «неопределённым поведением» в C? (Подсказка: что такое -(-128) в 8 битах?)
Исследовательский. Найдите максимальный int16. Какую физическую величину он ограничивал в Ariane 5?

Часть II. Точная арифметика за пределами фиксированного слова¶

§ 2. Арифметика произвольной точности: за пределами 64 бит¶

Провокация¶

В Python: 2 ** 1000 — число с 302 цифрами. Вычисляется мгновенно. В C: pow(2.0, 1024) уже переполняет double. А числа вроде 2**53 + 1 в float64 уже не различаются: соседние целые начинают сливаться.

Одна задача, разные ответы. Как Python хранит число с 302 цифрами в обычной оперативной памяти?

§ 2.1 Архитектура Python int: массив «цифр» в системе 2³⁰¶

Тип int в Python реализован не как фиксированная ячейка памяти, а как динамический массив. Python хранит большое целое число как последовательность 30-битных «цифр» (в документации CPython они называются digits) с основанием $2^{30}$ = 1 073 741 824. Это аналог столбчатой арифметики из начальной школы — только вместо «цифр» от 0 до 9 используются «цифры» от 0 до $2^{30}$-1, и столбцов может быть столько, сколько нужно.

Структура Python bigint¶

Число хранится как знак + массив d[0], d[1], ..., d[k-1], где 0 ≤ d[i] < $2^{30}$.

Значение числа: d[0] + d[1]·$2^{30}$ + d[2]·$2^{60}$ + ... + d[k-1]·$2^{30(k-1)}$

Число $2^{1000}$ хранится в ⌈1000/30⌉ = 34 элементах — всего ~136 байт памяти.

Сложение: O(n) операций. Умножение наивным алгоритмом: O(n²). Деление (алгоритм Кнута D): O(n²).

Анатолий Карацуба (1937–2008)¶

Советский и российский математик, специалист по теории чисел. Колмогоров выдвинул гипотезу: умножение двух n-значных чисел невозможно быстрее, чем за O(n²) операций, и организовал семинар для обсуждения. Карацуба, присутствовавший на семинаре, за одну неделю нашёл алгоритм с меньшей сложностью — опровергнув гипотезу. Работа была опубликована в 1962 году.

§ 2.2 Алгоритм Карацубы: «разделяй и властвуй» в арифметике¶

Наивное умножение двух n-значных чисел требует n² однозначных умножений. В 1960 году советский математик Анатолий Карацуба показал: это можно сделать быстрее.

Идея алгоритма Карацубы¶

Пусть a и b — числа с 2n цифрами. Представим: a = a₁·Bⁿ + a₀, b = b₁·Bⁿ + b₀, где B — основание.

Наивное произведение: a·b = a₁b₁·B²ⁿ + (a₁b₀ + a₀b₁)·Bⁿ + a₀b₀ — нужно 4 умножения.

Трюк Карацубы: p = a₁b₁, q = a₀b₀, r = (a₁+a₀)(b₁+b₀). Тогда a₁b₀ + a₀b₁ = r - p - q.

Итого: 3 умножения вместо 4. Плюс несколько сложений (они быстрые).

Сложность: T(n) = 3T(n/2) + O(n) → T(n) = O(n^log₂3) ≈ O(n¹·⁵⁸⁵). При n = $10^6$ это в ~1600 раз быстрее наивного.

§ 2.3 Сравнение алгоритмов умножения¶

Алгоритм	Сложность	Применение	Порог
Школьный	O(n²)	Python: < 70 цифр	< 70 цифр
Карацуба	O(n¹·⁵⁸⁵)	Python: 70–300 цифр	70–300 цифр
Тоом-Кук	O(n¹·⁴⁶⁵)	GMP: средние	> 300 цифр
Шёнхаге–Штрассен	O(n log n log log n)	GMP: миллионы цифр	> 10⁶ цифр
Harvey–Hoeven	O(n log n)	Теор. предел (2019)	N/A

Код: алгоритм Карацубы и быстрое возведение в степень¶

Что демонстрирует: рекурсивное умножение «разделяй и властвуй»; бинарное возведение в степень по модулю (основа RSA).

def karatsuba(a, b):
    if a < 1000 or b < 1000:
        return a * b
    n = max(len(str(a)), len(str(b)))
    m = n // 2
    B = 10 ** m
    a1, a0 = divmod(a, B)
    b1, b0 = divmod(b, B)
    p = karatsuba(a1, b1)
    q = karatsuba(a0, b0)
    r = karatsuba(a1 + a0, b1 + b0)
    return p * B**(2*m) + (r - p - q) * B**m + q
def fast_pow_mod(base, exp, mod):
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = result * base % mod
        base = base * base % mod
        exp >>= 1
    return result

Ваш браузер умножает 2048-битные числа прямо сейчас¶

Протокол HTTPS требует RSA-операций с числами длиной 617 десятичных цифр. Типичная операция: возвести 617-значное число в степень ~65537 по 617-значному модулю. Без быстрых алгоритмов (Карацуба + square-and-multiply) каждое соединение занимало бы секунды. С ними — единицы миллисекунд.

Задачи к §2¶

Базовый. Оцените: сколько байт занимает Python bigint, равный $2^{1000}$? Используйте sys.getsizeof().
Аналитический. Сравните скорость fast_pow_mod(2, 10**6, 10**9+7) и наивного метода. Измерьте через timeit.
Исследовательский. При каком размере чисел (в цифрах) алгоритм Карацубы начинает выигрывать у наивного? Проведите эксперимент.

§ 3. Дробные числа без потерь: рациональная арифметика¶

Провокация¶

Что такое 1/3 + 1/6? Любой школьник ответит: 1/2. Проверим:

Python: 1/3 + 1/6 → 0.49999999999999994

Почему не 0.5? Потому что 1/3 не имеет конечного двоичного представления. В памяти хранится не дробь 1/3, а «ближайшее представимое» число.

§ 3.1 Проблема: некоторые дроби не имеют конечного представления¶

Дробь p/q имеет конечное представление в системе с основанием b тогда и только тогда, когда все простые делители q делят b. В двоичной системе (b = 2): только дроби со знаменателем вида $2^k$ представляются точно.

Критерий точного двоичного представления¶

Дробь p/q (несократимая) имеет конечное двоичное представление тогда и только тогда, когда q = $2^k$ для некоторого натурального k.

Точные: 1/2 = 0.1₂, 1/4 = 0.01₂, 3/8 = 0.011₂

Неточные: 1/3, 1/5, 1/6, 1/7, 1/9, 1/10 = 0.0001100110011...₂

Следствие: числа 0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 0.8, 0.9 не представимы точно в double.

§ 3.2 Рациональная арифметика: хранить числитель и знаменатель¶

Решение: хранить рациональное число не как приближённое вещественное, а как пару целых (числитель, знаменатель). Все операции выполняются точно. Ошибок округления нет.

Рациональная арифметика: операции¶

Представление: q = p/r, где p ∈ ℤ, r ∈ ℕ, НОД(|p|, r) = 1.

Сложение: a/b + c/d = (a·d + b·c) / (b·d), затем сократить на НОД.

Умножение: a/b · c/d = (a·c) / (b·d), затем сократить на НОД.

Цена точности: числители и знаменатели растут. Сумма n дробей может дать числитель с тысячами цифр → неявно используется BigInt.

Код: Fraction и точная арифметика¶

Что демонстрирует: точное вычисление дробей без потери знаков; что именно хранится «под видом 0.1» в float64.

from fractions import Fraction
a = Fraction(1, 3)
b = Fraction(1, 6)
print(a + b)              # Fraction(1, 2) -- точно!
exact = sum(Fraction(1, n) for n in range(1, 101))
approx = sum(1/n for n in range(1, 101))
print(f'Ошибка float: {abs(float(exact) - approx):.2e}')
# Что реально хранится под видом '0.1'
print(Fraction(0.1))      # 3602879701896397/36028797018963968

3.3 Где рациональная арифметика применяется¶

Компьютерная алгебра (Mathematica, SymPy, Maple): символьные вычисления требуют точности. Финансы: Python Decimal, Java BigDecimal. Компиляторы: вычисление констант при компиляции. Верификаторы теорем (Coq, Lean, Isabelle). Компьютерная геометрия: предикаты «слева/справа от прямой» должны быть точными.

Patriot не перехватил Scud: 28 жертв из-за 0.1¶

Дахран, Саудовская Аравия, 25 февраля 1991 года. Иракская ракета Scud поразила казарму. 28 погибших, 98 раненых. Зенитный комплекс Patriot не произвёл пуск.

Причина: система отслеживала время в единицах 1/10 секунды. Это значение хранилось как 24-битное двоичное число. $1/10 = 0.0001100110011\ldots_2$ — бесконечная дробь. При обрезке до 24 бит: ошибка $\approx 9.54 \cdot 10^{-8}$ на каждые 0.1 с.

За 100 часов непрерывной работы накопилась ошибка: $0.0000000954 \cdot 100 \cdot 3600 \cdot 10 \approx 0.34$ секунды. Ракета Scud летит со скоростью ~1700 м/с. За 0.34 с — ~580 метров. Система искала цель в неверном секторе.

Если бы система отслеживала накопленную погрешность или периодически перезагружалась (как рекомендовала документация!) — катастрофы не было бы. Хранение времени в рациональной арифметике устранило бы ошибку полностью, но для встраиваемой системы реального времени 1991 года это было технически невозможно.

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

Задачи к §3¶

Базовый. Вычислите через Fraction точное значение 1 + 1/4 + 1/9 + ... + 1/n² для n = 100. Сравните с π²/6.
Аналитический. Запишите двоичное представление 0.1 с точностью 53 знаков (мантисса double).
Исследовательский. Если рациональная арифметика точная, почему её не используют везде? Подсказка: сравните sys.getsizeof() и скорость $10^6$ операций для float и Fraction.

Часть III. Float — механизм, не магия¶

§ 4. Вещественные числа: стандарт IEEE 754¶

Провокация¶

В Python: float('nan') == float('nan') даёт False. «Не число» не равно самому себе.

А float('inf') - float('inf') даёт nan, а не 0.

И -0.0 отличается от +0.0 на физическом уровне, хотя математически они равны.

Вещественная арифметика компьютера содержит объекты, которых нет в школьной математике. Разберём, откуда они берутся.

§ 4.1 Идея: число в форме ±m · 2ᵉ¶

Числа с плавающей точкой — это аналог научной нотации в двоичной системе. Так же, как физик пишет 3.14 $10^8$ м/с, компьютер хранит число в форме ±1.мантисса 2^порядок. Стандарт IEEE 754, принятый в 1985 году и обновлённый в 2008, определяет точный формат хранения для всех современных процессоров. Именно поэтому Python, C, Java, Rust — все дают одинаковый результат для float-вычислений.

Структура float32 (single precision, 32 бита)¶

Бит 31: знак (s). 0 = положительное, 1 = отрицательное.

Биты 30–23: порядок (e, 8 бит). Хранится со смещением 127: реальный порядок = e - 127.

Биты 22–0: мантисса (m, 23 бита). Хранит дробную часть числа 1.xxxxx. Целая «1» подразумевается.

Значение: $(-1)^s \cdot 1.m \cdot 2^{e-127}$. Диапазон: ±1.18·$10^{-38}$ ... ±3.4·$10^{38}$. Точность: ~7 десятичных знаков.

Структура float64 / double (64 бита)¶

Знак: 1 бит. Порядок: 11 бит, смещение 1023. Мантисса: 52 бита.

Значение: $(-1)^s \cdot 1.m \cdot 2^{e-1023}$. Диапазон: ±2.2·$10^{-308}$ ... ±1.8·$10^{308}$. Точность: ~15–16 знаков.

В стандартной реализации CPython все числа с плавающей точкой — это float64.

Разбор числа 0.375 = 3/8 в float32¶

Шаг 1. 0.375 > 0, знаковый бит s = 0.

Шаг 2. 0.375 = 0.011₂ = 1.1₂ $2^{-2}$, порядок e = -2. Хранится: e + 127 = 125 = 01111101₂.

Шаг 3. Мантисса: дробная часть 1.1₂ = 1₂ = 10000000000000000000000₂ (23 бита).

Итого: 0 01111101 10000000000000000000000₂

Проверка: $(-1)^0 \cdot 1.1_2 \cdot 2^{125-127} = 1.5 \cdot 2^{-2}$ = 0.375.

§ 4.2 Специальные значения: расширенная числовая система¶

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

Значение	Порядок e	Мантисса	Когда возникает
+Infinity	все 1	m=0, s=0	1.0/0.0, overflow
−Infinity	все 1	m=0, s=1	−1.0/0.0, −overflow
NaN	все 1	m≠0	0/0, sqrt(−1), inf−inf
−0.0	e=0	m=0, s=1	Предел слева
Денорм.	e=0	m≠0	Числа < 2⁻¹²⁶

IEEE 754 vs Python¶

Таблица выше описывает поведение на уровне стандарта IEEE 754 (и, соответственно, C, NumPy, Rust). В Python на уровне языка поведение отличается: 1.0 / 0.0 вызывает ZeroDivisionError (а не даёт +Infinity), math.sqrt(-1) вызывает ValueError (а не даёт NaN). Через NumPy (np.float64(1.0) / np.float64(0.0)) вы увидите именно IEEE-поведение.

Поведение специальных значений¶

На уровне IEEE/NumPy: inf + 1 → inf, inf - inf → nan, inf · 0 → nan

В Python math: 1.0 / 0.0 → ZeroDivisionError (не +Infinity!). Через NumPy: np.float64(1.0) / np.float64(0.0) → inf.

nan == nan → False — NaN не равен ничему, включая себя. Это поведение одинаково и в Python, и в C, и в NumPy.

-0.0 == 0.0 → True (равны при сравнении), но 1/(-0.0) → -inf через NumPy (в Python math — ZeroDivisionError). На уровне битов -0.0 и +0.0 различимы.

§ 4.3 Машинное эпсилон: главная характеристика точности¶

Машинное эпсилон eps — наименьшее положительное число такое, что 1.0 + eps ≠ 1.0.

float32: eps ≈ $1.19 \cdot 10^{-7}$. Точность: ~7 десятичных знаков.

float64: eps ≈ $2.22 \cdot 10^{-16}$. Точность: ~15–16 знаков.

Физический смысл: две соседние float-точки вблизи 1.0 отличаются на eps. Любое вещественное число x округляется до ближайшего float с относительной ошибкой ≤ eps/2.

Python: import sys; sys.float_info.epsilon → 2.220446049250313e-16

§ 4.4 Неравномерное распределение float на числовой оси¶

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

Float-точки не равноудалены¶

Между 1.0 и 2.0 в float32 ровно $2^{23}$ = 8 388 608 представимых чисел. Шаг $\approx 1.19 \cdot 10^{-7}$.

Между 2.0 и 4.0 — тоже $2^{23}$, но шаг вдвое больше: $\approx 2.38 \cdot 10^{-7}$.

Между 1024 и 2048 — снова $2^{23}$, но шаг $= 1/8192 \approx 1.2 \cdot 10^{-4}$. В 1000 раз хуже, чем вблизи 1.

Между $10^{15}$ и $10^{15} + 1$ в float64 нет ни одного представимого числа. То есть 10**15 + 0.5 == 10**15 в double-арифметике.

Что отличает 0.1 + 0.2 от 0.3 — теперь по битам¶

Возвращаемся к первому вопросу из начала главы. Тождество 0.1 + 0.2 == 0.3 ложно в любом языке с IEEE 754. Причина теперь видна на уровне битов: ни 0.1, ни 0.2 не имеют конечного двоичного представления (их знаменатель — не степень двойки). Оба числа округляются до ближайших float64-точек; сумма этих приближений отличается от приближения суммы ровно на 1 ULP — единицу в последнем разряде мантиссы. Это и есть тот бит, о котором спрашивал первый вопрос. Инженерное следствие просто и обязательно: сравнение float-результатов — всегда через допуск; math.isclose(a, b) либо abs(a - b) <= atol + rtol * scale, и никогда не ==.

Код: побитовое устройство float64 и машинное эпсилон¶

Что демонстрирует: разложение float64 на знак, порядок и мантиссу; экспериментальное нахождение eps.

import struct, sys
def float_to_bits(x):
    bits = struct.pack('>d', x)
    n = int.from_bytes(bits, 'big')
    s = (n >> 63) & 1
    e = (n >> 52) & 0x7FF
    m = n & ((1 << 52) - 1)
    return f'sign={s} exp={e}(={e-1023}) mantissa={m:052b}'
print(float_to_bits(0.1))
print(float_to_bits(0.3))
print(float_to_bits(0.1 + 0.2))  # отличается от 0.3!
# Найти машинное эпсилон экспериментально
eps = 1.0
while 1.0 + eps != 1.0:
    eps /= 2
eps *= 2
print(f'eps = {eps:.2e}')        # ~2.22e-16

Задачи к §4¶

Базовый. Найдите машинное эпсилон для float32 и float64 экспериментально. Сравните с теоретическим.
Аналитический. Объясните: почему (1.0 + 1e-16) == 1.0 в Python? Вычислите: при каком x ещё (1.0 + x) != 1.0?
Исследовательский. Найдите наименьшее целое n такое, что float(n) == float(n+1). Объясните через структуру float64.

§ 5. Как накапливается ошибка: три механизма¶

Провокация¶

sum(0.1 for _ in range(10)) → 0.9999999999999999 (не 1.0!)

Сложили ровно 10 раз по 0.1. Откуда ошибка в последнем знаке? Это не единственный сценарий. Бывает, что ошибка не в 16-м знаке, а в первом.

Механизм 1: ошибка представления (roundoff error)¶

Правило: $\mathrm{fl}(x) = x \cdot (1 + \delta)$, где $|\delta| \le \text{eps}/2$.

Для float64: $|\delta| \le 1.11 \cdot 10^{-16}$. Мало.

Для $n$ операций: ошибка $O(n \cdot \text{eps})$. При $n = 10^6$ — $\sim 10^{-10}$. Терпимо.

Проблема начинается, когда другие механизмы усиливают эту ошибку.

Механизм 2: катастрофическое вычитание (cancellation)¶

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

Пусть $a \approx b$. Оба известны с относительной погрешностью eps.

$a - b =$ маленькое число. Абсолютная погрешность: $\sim \text{eps} \cdot |a|$.

Относительная погрешность результата: $\text{eps} \cdot |a| / |a - b|$.

Если $|a - b| \ll |a|$, то относительная погрешность $\gg 1$. Все значащие цифры потеряны.

Пример: $a = 1.000000001$, $b = 1.000000000$ (с точностью $\pm 10^{-16}$). Разность $= 10^{-9}$. Относительная ошибка $\sim 10^{-7}$ — потеряно 7 знаков!

Катастрофическое вычитание в квадратном уравнении¶

Задача: корни $x^2 - 10000x + 1 = 0$. Дискриминант $D \approx 10^8$, $\sqrt{D} \approx 9999.9998$.

$x_1 = (10000 + 9999.9998)/2 \approx 9999.9999$ ← вычисляется хорошо.

$x_2 = (10000 - 9999.9998)/2 \approx 0.0001$ ← катастрофическое вычитание!

Исправление: $x_1 = (-b - \mathrm{sign}(b) \cdot \sqrt{D}) / 2a$, $x_2 = c / (a \cdot x_1)$ — теорема Виета устраняет вычитание.

Механизм 3: число обусловленности¶

Некоторые задачи чувствительны к малым изменениям входных данных. Такие задачи называются плохо обусловленными. Компьютер не виноват — математика такова.

$K =$ отношение относительного изменения результата к относительному изменению входных данных.

$K \sim 1$: хорошо обусловленная задача. $K \gg 1$: плохо обусловленная.

Для системы $Ax = b$: $K(A) = \|A\| \cdot \|A^{-1}\|$. Правило: $K = 10^k \Rightarrow$ теряем $\sim k$ знаков.

Важно: плохая обусловленность — свойство задачи, не алгоритма. Даже идеальный алгоритм не поможет.

Полином Уилкинсона¶

$p(x) = (x-1)(x-2)\ldots(x-20)$. Корни: $1, 2, \ldots, 20$.

Уилкинсон изменил коэффициент при $x^{19}$ на $-210 - 2^{-23}$ (изменение в 8-м знаке!).

Результат: корни разлетелись в комплексную плоскость. 10 комплексных пар вместо 20 вещественных. Число обусловленности $\sim 10^{13}$.

Вывод: задача нахождения корней полинома по коэффициентам — принципиально плохо обусловленная.

Код: cancellation, устойчивое квадратное уравнение, обусловленность¶

Что демонстрирует: потеря знаков при вычитании близких чисел; исправление через ряд Тейлора и теорему Виета; рост cond(H) для матриц Гильберта.

import numpy as np
# Катастрофическое вычитание
x = 1e-8
direct = x - np.sin(x)            # плохой способ
taylor = x**3/6 - x**5/120        # хороший способ
print(f'Прямое: {direct:.6e}, Тейлор: {taylor:.6e}')
# Устойчивое квадратное уравнение
def stable_quadratic(a, b, c):
    D = b**2 - 4*a*c
    if D < 0: return None
    sqrt_D = np.sqrt(D)
    if b >= 0: x1 = (-b - sqrt_D) / (2*a)
    else:      x1 = (-b + sqrt_D) / (2*a)
    x2 = c / (a * x1)
    return x1, x2
# Число обусловленности матрицы Гильберта
H = np.array([[1.0/(i+j-1) for j in range(1,11)] for i in range(1,11)])
print(f'cond(H10) = {np.linalg.cond(H):.2e}')  # ~10^13

Задачи к §5¶

Базовый. Вычислите x - sin(x) при x = $10^{-8}$ двумя способами. Используйте mpmath для эталона. Сколько верных знаков?
Аналитический. Сравните stable_quadratic и наивную формулу для x² - $10^8$x + 1 = 0.
Исследовательский. Вычислите cond(H(n)) для n = 5, 8, 10, 12. Что происходит?

Часть IV. Численная инженерия: как не сделать ошибку¶

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

§ 6. Алгоритмы и точность: устойчивые против неустойчивых¶

Провокация¶

Два алгоритма вычисляют одно: сумму n чисел. Первый — наивный. Второй — алгоритм Кэхэна: четыре строки дополнительного кода.

Математически они тождественны. Численно — разница в точности в миллион раз при n = $10^6$.

§ 6.1 Численная устойчивость: определение¶

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

Численно неустойчивый: ошибка может расти экспоненциально — результат теряет все значащие цифры.

Прямая устойчивость (forward stability): ошибка результата мала — алгоритм возвращает ответ, близкий к точному.

Обратная устойчивость (backward stability): результат является точным решением слегка возмущённой задачи. Более сильное и полезное свойство.

Ключевой факт: плохая обусловленность + устойчивый алгоритм = приемлемый результат. Хорошая обусловленность + неустойчивый алгоритм = потеря точности.

Уильям Кэхэн (род. 1933)¶

Канадский математик, профессор UC Berkeley. Главный архитектор стандарта IEEE 754. Лауреат премии Тьюринга (1989) — за работу по стандарту IEEE 754 и численным методам. Алгоритм компенсированного суммирования — 1965 год.

§ 6.2 Алгоритм Кэхэна: компенсированное суммирование¶

Наивное суммирование: при прибавлении маленького слагаемого к большой сумме «маленькая» часть теряется при округлении. Кэхэн эту потерю отслеживает и компенсирует.

Наивное суммирование n чисел: ошибка O(n · eps · max|xᵢ|). При n = $10^6$ теряем 6 знаков.

Алгоритм Кэхэна: ошибка O(eps · |результат|) независимо от n. Практически машинная точность.

Идея: на каждом шаге вычисляем «потерянную» часть и передаём на следующий шаг как компенсацию.

Код: алгоритм Кэхэна и сравнение с наивным¶

Что демонстрирует: компенсация потери при суммировании; разница в точности на $10^6$ слагаемых.

def kahan_sum(numbers):
    total = 0.0
    compensation = 0.0
    for x in numbers:
        y = x - compensation
        t = total + y
        compensation = (t - total) - y  # что потеряно
        total = t
    return total
import math
from fractions import Fraction
n = 1_000_000
data = [0.1] * n
exact = float(Fraction(1, 10) * n)
print(f'Наивное:  ошибка {abs(sum(data)-exact):.2e}')
print(f'Кэхэн:   ошибка {abs(kahan_sum(data)-exact):.2e}')
print(f'fsum:     ошибка {abs(math.fsum(data)-exact):.2e}')

§ 6.3 Решение систем линейных уравнений: ведущий элемент¶

Метод Гаусса — главный алгоритм линейной алгебры. Но в наивной реализации скрыта численная проблема: деление на малый диагональный элемент усиливает ошибки.

Метод Гаусса: устойчивый против неустойчивого¶

Без выбора ведущего элемента: O(n³). Неустойчив. В продакшн-коде не использовать.

С частичным выбором ведущего (partial pivoting): O(n³). Устойчив для большинства задач. Стандарт: LAPACK, scipy.linalg.solve.

Для более сложных случаев существуют LU-разложение, QR-разложение и SVD (сингулярное разложение) — каждый со своей областью применения.

Sleipner A: 700 миллионов долларов и 47% ошибки¶

23 августа 1991 года. Норвегия. Бетонное основание нефтяной платформы Sleipner A разрушилось при контролируемом погружении на рабочую глубину в Северном море. Ущерб: около 700 миллионов долларов.

Причина: численная ошибка в конечно-элементном расчёте прочности бетонных перегородок. Слишком грубая сетка усилила ошибку в оценке сдвиговых напряжений. Расчётная нагрузка была занижена примерно на 47%. Перегородки треснули под давлением воды.

Урок прост: верификация численного алгоритма — не академический вопрос, а вопрос безопасности. То, что в учебнике выглядит как «выбор разрешения сетки», в инженерной практике может стоить сотен миллионов долларов и человеческих жизней.

Задачи к §6¶

Базовый. Сравните наивное суммирование и Кэхэна для n = $10^6$ случайных float64 из [0,1]. Fraction как эталон.
Аналитический. Создайте матрицу Гильберта H(n), n = 5, 8, 10, 12. Решите Hx = b (x = вектор единиц). График ошибки от n.
Исследовательский. Объясните: что такое «невязка» и почему она лучше характеризует качество решения, чем ‖x_комп - x_точн‖?

Углубление¶

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

§ 7. Как оценивать точность вычислений¶

Провокация¶

Нельзя писать: if computed_value == exact_value: — это почти никогда не сработает для float. Даже if result == 0.0 ненадёжно.

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

§ 7.1 Абсолютная и относительная погрешность¶

Абсолютная: Δx = |x_computed - x_exact|. Размерная. Зависит от масштаба.

Относительная: δx = |x_computed - x_exact| / |x_exact|. Безразмерная.

Когда абсолютную: задача имеет физический масштаб. Когда относительную: сравнение чисел разного масштаба.

Проблема: относительная не определена при x_exact = 0. Используйте смешанную: max(|x|, |y|) в знаменателе.

§ 7.2 ULP — единицы на последнем месте¶

При первом чтении достаточно вынести § 7.1, § 7.3 и инженерный чеклист (§ 7.5). Разделы 7.2 и 7.4 — более тонкие инструменты для углублённого анализа.

ULP (Unit in the Last Place) — расстояние от числа $x$ до следующего представимого float. Два числа различаются на 1 ULP, если между ними нет ни одного float.

math.ulp(x) в Python 3.9+: расстояние от x до следующего float.

Результат «правильно округлён», если отличается от точного менее чем на 0.5 ULP.

math.ulp(1.0) ≈ $2.22 \cdot 10^{-16}$ — это в точности машинное эпсилон. Почему? По определению eps — расстояние от 1.0 до следующего float, то есть ULP(1.0).

§ 7.3 Правильное сравнение float¶

Для результатов вещественных вычислений: почти никогда не сравнивайте a == b напрямую; используйте допуск.

Абсолютное: abs(a - b) < atol — работает при числах близких к нулю.

Относительное: abs(a - b) / max(abs(a), abs(b)) < rtol.

Комбинированное (numpy): abs(a - b) ≤ atol + rtol · max(abs(a), abs(b)).

Рекомендация: rtol = $10 \cdot \text{eps} \approx 10^{-15}$ для одиночных операций, rtol = $10^{-10}$ для длинных вычислений.

§ 7.4 Интервальная арифметика¶

Радикальный подход: вместо числа x хранить интервал [x_min, x_max], гарантированно содержащий точное значение. Каждая операция обновляет интервал.

[a₁,a₂] + [b₁,b₂] = [a₁+b₁, a₂+b₂].

[a₁,a₂] [b₁,b₂] = [min(...), max(...)] по всем парам.

Применение: доказательство гипотезы Кеплера об оптимальной упаковке шаров (Томас Хейлз, 2005) использовало интервальную арифметику для компьютерной верификации. Также: управление беспилотниками, финансовые расчёты.

Цена: замедление в 10–100 раз. Ширина интервала может расти быстро.

Код: сравнение float и ULP-расстояние¶

Что демонстрирует: стандартный практический способ сравнения float через math.isclose(); демонстрационный подсчёт ULP-расстояния для положительных конечных float (не универсальная реализация).

import math, struct
# Стандартный способ: math.isclose()
# Реализует комбинированный критерий: abs(a-b) <= atol + rtol * max(abs(a), abs(b))
print(math.isclose(0.1 + 0.2, 0.3))              # True (rel_tol=1e-9)
print(math.isclose(1e-15, 2e-15, abs_tol=1e-14))  # True (вблизи нуля нужен abs_tol)
# ULP-расстояние (демонстрация для положительных конечных float)
def ulp_distance(a, b):
    assert a > 0 and b > 0, "только положительные конечные"
    def to_int(x):
        n = int.from_bytes(struct.pack('>d', x), 'big')
        return n
    return abs(to_int(a) - to_int(b))
print(ulp_distance(0.1 + 0.2, 0.3))  # несколько ULP

§ 7.5 Инженерный чеклист¶

Ситуация	Действие
Начинаете вычисления?	Оцените cond. K > 10¹² → спец. методы
Суммируете много чисел?	Кэхэн или `math.fsum()`. Не `sum()`.
Сравниваете вычисленные вещественные результаты?	Почти никогда не `==`. `math.isclose()` или комбинированный допуск (`atol`, `rtol`).
Вычитаете близкие числа?	Переформулируйте алгебраически.
Решаете СЛАУ?	С pivoting. Проверьте cond(A).
Нужна гарантия?	`Fraction` или `mpmath`.
Физическая задача?	Точность ≥ точности измерений.

Теорема Годунова: математический предел точности¶

В 1961 году Сергей Годунов доказал: не существует линейного монотонного численного метода для гиперболических уравнений точнее первого порядка. Любой «более точный» линейный метод порождает осцилляции вблизи разрывов (ударные волны, фронты горения).

Выход: нелинейные ограничители потока, адаптивные сетки, WENO-схемы.

Урок: иногда «более точный алгоритм» даёт физически неверный результат, а «грубый» — правильный. Точность и устойчивость — разные качества.

Задачи к §7¶

Базовый. Протестируйте math.isclose(a, b) на парах: (1.0, 1.0+1e-10), (0.0, 1e-300), (1e300, 1e300+1e284). Для каких пар нужен abs_tol, а для каких достаточно rel_tol?
Аналитический. Используя mpmath (50 знаков), вычислите Σ 1/n² для n = 1..10000. Сравните с π²/6. Сколько знаков верно в float64?
Исследовательский. (Проект) Решение квадратного уравнения с явным анализом ошибок: для каждого корня — оценка погрешности через cond.

§ 8. Дерево решений: какой числовой тип выбрать?¶

Это кульминация главы. Все предыдущие параграфы вели к одному практическому вопросу: имея конкретную задачу, какой числовой тип взять?

Decision tree: выбор числового типа¶

1. Число целое?

→ Диапазон заранее известен и $\le 64$ бит? → int32 / int64.

→ Диапазон непредсказуем или нужна криптография? → Python int (BigInt).

2. Число дробное и нужна точность?

→ Финансы, символьные вычисления, предикаты геометрии? → Fraction.

→ Финансы с фиксированным числом десятичных знаков? → Decimal.

3. Число дробное и нужна скорость?

→ GPU, графика, нейросети, потоковые данные? → float32.

→ Наука, инженерия, общие вычисления? → float64.

4. Нужна гарантированная точность или верификация? → Интервальная арифметика.

Правило: начинайте с самого простого типа, который решает задачу. Переходите к более тяжёлому только когда есть конкретная причина — не «на всякий случай».

Сводная таблица: числа в компьютере¶

Тип	Размер	Диапазон	Точность	Применение	Катастрофа
int8/uint8	1 Б	−128..127 / 0..255	точно	Байты, пиксели	—
int16/uint16	2 Б	−32768..32767	точно	Порты, аудио	Ariane 5*
int32/uint32	4 Б	±2.1·10⁹	точно	Стандартный int	—
int64/uint64	8 Б	±9.2·10¹⁸	точно	long, Unix-время	—
Python int	перем.	неограничен	точно	BigInt, крипто	—
float32	4 Б	±3.4·10³⁸	~7 знаков	GPU, графика	—
float64	8 Б	±1.8·10³⁰⁸	~16 знаков	Наука, инженерия	Patriot
Fraction	перем.	точно	точно	Символьн. выч.	—
Decimal	перем.	настр.	настр.	Финансы	—
mpmath	настр.	произвольно	произвольно	Верификация	—
Ariane 5: ошибка при преобразовании float64 → int16.

§ 9. Проектная задача: спасите вычисление¶

Задание¶

Дано выражение для вычисления:

\[f(x) = \frac{1 - \cos x}{x^2}, \quad x \to 0.\]

Точное значение предела: $1/2$.

import math
from mpmath import mp, mpf, cos as mpcos
mp.dps = 50  # 50 знаков для эталона
def f_direct(x):
    return (1 - math.cos(x)) / x**2
def f_trig(x):
    half = x / 2
    sinc = math.sin(half) / half if half != 0 else 1.0
    return sinc**2 / 2
def f_taylor(x):
    return 0.5 - x**2/24 + x**4/720
def f_exact(x):
    xm = mpf(x)
    return float((1 - mpcos(xm)) / xm**2)
for e in [1, 4, 8, 12]:
    x = 10.0 ** (-e)
    exact = f_exact(x)
    print(f'x=1e-{e}: direct err={abs(f_direct(x)-exact)/exact:.2e}, '
          f'trig err={abs(f_trig(x)-exact)/exact:.2e}, '
          f'taylor err={abs(f_taylor(x)-exact)/exact:.2e}')

Проект: катастрофическое cancellation и его лечение¶

Часть 1. Вычислите f(x) напрямую для $x = 10^{-1}, 10^{-4}, 10^{-8}, 10^{-12}$. Сравните с $0.5$. Постройте таблицу относительных ошибок. Что происходит?
Часть 2. Объясните механизм: при малых $x$ числитель $(1 - \cos x)$ вычисляет разность двух чисел, близких к 1. Это катастрофическое вычитание.
Часть 3. Перепишите формулу так, чтобы уменьшить ошибку в 1000 раз и более. Два способа:

Способ A (алгебраический): используйте тождество $1 - \cos x = 2 \sin^2(x/2)$, тогда $f(x) = 2 \sin^2(x/2) / x^2 = (\sin(x/2) / (x/2))^2 / 2$. При малых $x$ $\sin(x/2)/(x/2) \to 1$, и cancellation исчезает.

Способ B (ряд Тейлора): $\cos x \approx 1 - x^2/2 + x^4/24$, тогда $(1 - \cos x)/x^2 \approx 1/2 - x^2/24$. Вычитание устранено.

Часть 4. Постройте таблицу ошибок для обоих способов. При каком $x$ прямая формула теряет все знаки, а исправленная — нет?
Часть 5. Напишите функцию safe_f(x), которая автоматически выбирает способ вычисления в зависимости от $|x|$.

Итог: от транзистора до гарантированной точности¶

Эта глава проследила путь от транзистора до интервальной арифметики. От «сколько бит» до «как гарантировать правильность результата».

Ключевой сдвиг мышления: компьютер не вычисляет «точно». Он вычисляет «достаточно точно при правильных условиях». Понять эти условия — значит стать инженером, а не просто пользователем калькулятора.

Главный вывод главы¶

Точность вычислений — это инженерный компромисс. Целые числа: точны, но ограничены диапазоном. BigInt: точны, но медленнее. Рациональные: точны, но растут бесконтрольно. Float: быстрые и компактные, но приближённые. Понимать этот компромисс — знать, какой инструмент взять для каждой задачи.

Ответы на три вопроса из начала главы¶

0.1 + 0.2 == 0.3 → False. Числа 0.1 и 0.2 не имеют конечного двоичного представления (знаменатель не степень двойки). Сумма приближений отличается от приближения суммы на 1 ULP — конкретный бит мантиссы, который мы увидели через float_to_bits().
int в Python — BigInt: динамический массив 30-битных «цифр», никогда не переполняется. int в Java — int32: фиксированная ячейка, арифметика по модулю $2^{32}$, гарантированное переполнение. int в C (signed) — фиксированная ячейка, но переполнение формально не определено стандартом (undefined behavior); unsigned — по модулю $2^n$.
Ariane 5: преобразование float64 → int16 без проверки диапазона. Горизонтальная скорость превысила 32 767 — максимум int16. Необработанное исключение → диагностические данные вместо полётных → ракета разрушилась.

Мост к следующей главе¶

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

Этим мы займёмся в главе «Множества, логика и язык математики».

Что за пределами этой главы¶

Числа с фиксированной точкой (fixed-point): компромисс между integer и float, используется в DSP и встраиваемых системах.

SIMD и векторные вычисления: параллельная float-арифметика на уровне CPU и GPU.

Символьные вычисления: SymPy, Mathematica — работа с формулами, не числами.

Верификация численных программ: Frama-C, Why3 — формальное доказательство корректности float-программ.

Квантовые вычисления: иная модель числовой точности — вероятностные результаты.