Длинная арифметика

Факториал тысячи — число из 2568 цифр. В 64-битное машинное слово помещается девятнадцать, а Python печатает 1000! целиком и мгновенно. Внутри он хранит его как 285 «цифр» по 30 бит и умножает их так же, как мы умножаем столбиком, только каждая «цифра» — от нуля до 1 073 741 823. Ниже n! растёт в такой памяти умножение за умножением: каждая клетка — одна цифра, цвет — её значение. Переключите основание, а потом нажмите «Напечатать»: окажется, что вывести число на экран бывает дороже, чем его посчитать.

Чёрные клетки в начале — нули на конце числа: у 1000! по основанию 2³⁰ их 33, по основанию 10⁹ — 27, а в десятичной записи нулей 249. Посмотрите на печать. По основанию 10⁹ каждая клетка — готовые девять десятичных цифр, и печать — один проход по памяти. По основанию 2³⁰ десятичных цифр в памяти нет вовсе: чтобы получить последние девять, нужно разделить всё число на 10⁹, потом частное — ещё раз, и так 286 раз. Каждое деление проходит по всем клеткам, и для 285 клеток выходит больше сорока тысяч делений клеток — квадрат длины.

Идея

Длинное число — массив «цифр» по основанию B, младшая первой: x = d₀ + d₁B + d₂B² + … Сложение, вычитание, умножение и деление — школьные алгоритмы в столбик, только цифра теперь — машинное число. Основание берут как можно больше, но так, чтобы промежуточные результаты — сумма двух цифр с переносом и произведение двух цифр с переносом — помещались в машинное слово двойной ширины.

Сложение — проход от младшей цифры к старшей: s = a_i + b_i + перенос, цифра результата — s mod B, новый перенос — s div B, то есть 0 или 1. Умножение на короткое число — тоже один проход, перенос меньше множителя. Умножение двух длинных чисел столбиком — каждая цифра на каждую, O(n·m) умножений цифр. Насколько велико может быть основание, решает одно неравенство.

Отсюда выбор оснований. Олимпиадное B = 10⁹: B² = 10¹⁸ помещается в знаковые 64 бита, а 10¹⁰ уже нет — 10²⁰ больше 2⁶⁴. И печатать такое число легко: каждая цифра — девять десятичных знаков. Двоичные основания экономнее: при 2³² в деле каждый бит слова, произведение двух цифр помещается в 64 бита, сдвиги и побитовые операции почти бесплатны. Go, GMP и BigInt в JavaScript идут дальше и отдают цифре всё 64-битное слово: произведение двух таких цифр занимает 128 бит, но процессор и так выдаёт их одной командой — старшую и младшую половины, в Go это bits.Mul64.

Python выбрал 30 бит. В комментарии к исходному коду CPython перечислены условия: в типе цифры должно помещаться 2B − 1 — сумма двух цифр без переполнения, в типе двойной цифры — B² − 1, а формат сохранения marshal требует, чтобы длина цифры была кратна 15. Подходят 15 и 30 бит; на 64-битных машинах цифра — 30 бит в 32-битной ячейке, произведение — в 64-битной.

Код

Длинные числа по основанию 10⁹, как их пишут на олимпиадах: разбор строки, печать, сложение с подсчётом переносов, умножение на короткое число, столбик и деление на короткое. Программа считает 1000! и 2¹⁰⁰⁰ и сверяет всё со встроенной арифметикой языка.

1000! has 2568 digits in 286 base-10^9 digits, 249 zeros at the end
it starts with 40238726007709377354 and its last nonzero digits are 0027753472
2^1000 has 302 digits; 8063 carries in 17212 digit additions
(10^1000 - 1) * 3^2000: 112 x 107 digits -> 218, mod 1000000007 = 817895449

При удвоениях перенос случается в 47 % сложений цифр — когда цифра не меньше половины основания. Длинные цепочки переносов, как в 999…9 + 1, бывают, но редки, и сложение длинных чисел стоит O(n). Деление на короткое число идёт в обратную сторону, от старшей цифры: остаток умножается на B и прибавляется к следующей цифре, как при делении уголком.

Как это устроено в Python и Go

Целое число в Python — всегда длинное. Заголовок объекта занимает 24 байта, в последнем его поле — число цифр и знак, дальше идут сами цифры по 4 байта. Поэтому sys.getsizeof для числа из k цифр — 24 + 4k: у 1000! это 285 цифр и 1164 байта, ровно столько показывает первая сцена. Параметры видны прямо из языка: sys.int_info сообщает bits_per_digit=30 и sizeof_digit=4. Маленькие числа — тоже длинные, просто с одной цифрой, а сложение двух таких Python проверяет отдельной быстрой веткой.

В Go длинное число — big.Int: знак и срез машинных слов nat, младшее слово первым. Внутренние циклы — «сложить два вектора слов», «умножить вектор на слово и прибавить к другому» — написаны на ассемблере под каждый процессор, а bits.Mul64 и bits.Add64 компилируются в одну-две команды. Поэтому переписывать math/big вручную незачем, а на длинных числах она сама переключается со столбика на умножение Карацубы.

Деление: угадать цифру частного

Деление длинного на длинное — тот же уголок, но цифру частного нельзя подобрать перебором: их B, около миллиарда. Её угадывают по старшим цифрам: делят две старшие цифры остатка на старшую цифру делителя. Догадка может оказаться больше настоящей цифры, и тогда её уменьшают, проверяя по следующей цифре. Дональд Кнут в алгоритме D показал, что если заранее домножить оба числа так, чтобы старшая цифра делителя была не меньше B/2, догадка превосходит настоящую цифру частного не больше чем на 2. Доказательство и все детали — в «Искусстве программирования», том 2, § 4.3.1; здесь важен итог: деление столбиком, как и умножение, стоит O(n·m) операций с цифрами.

Для очень длинных чисел деление сводят к умножению. Обратную величину 1/b находят методом Ньютона — каждая итерация удваивает число верных цифр, — а частное получают умножением. Время деления тогда того же порядка, что и умножения. Другой путь — рекурсивное деление половинами, по той же схеме «разделяй и властвуй», что у Карацубы; так делят длинные числа GMP и math/big.

Печать: где прячется квадрат

Двоичное основание удобно для арифметики и неудобно для людей. Чтобы напечатать число, его переводят в десятичную запись, и наивный способ — тот, что в первой сцене: отщипывать снизу по девять цифр делением на 10⁹. Каждое деление проходит по всему числу, делений — по одному на девять цифр, итого O(n²). На числе в сто тысяч цифр это уже четверть секунды на Python, на миллионе — полминуты.

Быстрее — «разделяй и властвуй». Делим число на 10^9·2k — степень, близкую к корню из числа, — и получаем две половины по числу цифр: старшую и младшую. Каждую переводим тем же способом, младшую дополняем ведущими нулями. Степени 10⁹, 10¹⁸, 10³⁶, … считаются возведением в квадрат и хранятся в таблице. Деление половин стоит столько же, сколько умножение такой длины, M(n), уровней рекурсии — log n, и всё вместе — O(M(n) log n) по основной теореме.

Замер: число из d цифр переводится в строку тремя способами — нашим наивным, нашим «пополам» и встроенным, — и строка разбирается обратно в число.

Наивный способ растёт с наклоном 2: в десять раз длиннее — в сто раз дольше, на ста тысячах цифр 0,35 секунды на Python и 0,29 на Go. «Пополам» на тех же ста тысячах в 4 раза быстрее на Python и в 15 раз на Go, и разрыв растёт. Встроенные способы ведут себя по-разному. В Go String() устроен как «пополам» — в исходном коде math/big записано, что рекурсия ускоряет перевод вчетверо на средних длинах и примерно в 30 раз на двадцати тысячах цифр, — а вот разбор строки SetString почти квадратичный: на трёхстах тысячах цифр он уже втрое медленнее печати. В Python и печать, и разбор растут медленнее квадрата: начиная с версии 3.12 очень длинные числа CPython переводит алгоритмами «разделяй и властвуй» из модуля _pylong.

Сколько стоят операции

Памяти — O(n) на число, и её легко оценить заранее: n! содержит около log₁₀ n! ≈ n(log₁₀ n − 0,434) десятичных цифр по формуле Стирлинга, а 30-битных цифр Python — в log₁₀ 2³⁰ ≈ 9,03 раза меньше. А вот что умеет движок вашего браузера — BigInt устроен в нём по-своему. Нажмите «Замерить»: пять операций на числах от тысячи до трёхсот тысяч цифр.

Сложение растёт линейно, остальные операции — медленнее квадрата: в Chrome умножение у нас вышло с наклоном около 1,2, деление и печать — около 1,4–1,5. Разные движки JavaScript здесь заметно расходятся, поэтому вывод один: если программа работает с числами в сотни тысяч цифр, нужно мерить на той платформе, где она будет работать.

Где это работает

• Python, Go, JavaScript, Java. В Python каждое целое — длинное; в Go — math/big; в JavaScript — тип BigInt, появившийся в стандарте 2020 года; в Java — BigInteger. На олимпиадах на C++ длинную арифметику по основанию 10⁹ по-прежнему пишут сами.
• Криптография. Ключ RSA на 2048 бит — это 69 цифр Python или 32 слова Go, и всё шифрование — возведение в степень по модулю таких чисел. Для неё библиотеки вроде OpenSSL и GMP держат отдельные быстрые циклы умножения по модулю.
• Рекордные числа. В октябре 2024 года проект GIMPS нашёл простое число 2^{136 279 841} − 1 из 41 024 320 цифр. Проверка такого числа — 136 миллионов возведений в квадрат по модулю числа из 136 миллионов бит. Столбиком каждое стоило бы триллионы умножений слов, и проверка растянулась бы на тысячи лет; с умножением через быстрое преобразование Фурье она становится посильной. На этот раз считали тысячи графических процессоров.

История

Школьные алгоритмы с длинными числами старше компьютеров, а для машин их систематически описал Дональд Кнут во втором томе «Искусства программирования» в 1969 году: алгоритмы A, S, M и D — сложение, вычитание, умножение и деление — с доказательствами и анализом. Библиотеку GMP Торбьёрн Гранлунд начал в январе 1991 года и выпустил в сентябре; сегодня на ней держатся системы компьютерной алгебры и многие рекордные вычисления. В Python длинные целые были с первых версий, а с Python 3.0 (2008) остался единственный тип int. В Go пакет math/big вошёл в первую стабильную версию языка в 2012 году.

Задачи

В первой задаче нужно сосчитать переносы при сложении двух чисел по миллиону цифр — среди тестов первый миллион цифр π и e. Во второй — сколько нулевых младших ячеек у n! по основанию 2³⁰, 10⁹ или любому другому, когда сам факториал не поместится ни в какую память.