Нахождение первых 1000 знаков после запятой в записи числа e

Полагаю, читателю не нужно объяснять, что такое число e и какую роль в математике оно играет. Напомню лишь, что оно является иррациональным. Это означает, что его можно представить в виде бесконечной десятичной дроби, причём, эта дробь является непериодической.

Достаточно давно в одном из задачников по программированию я встретил 2 задачи, в одной из которых требовалось найти 100 знаков после запятой в записи числа e, а в другой — числа π. Тогда я задумался о том, как можно решить эти задачи, но ничего путного не придумал. А сейчас решил вернуться к первой из них (возможно, потом вернусь и ко второй). Но, благодаря статье из Википедии, о которой мы ещё поговорим, решил не ограничиваться 100 знаками, а найти сразу 1000.

Давайте теперь сформулируем саму задачу.

Найти первые 1000 знаков после запятой в десятичной записи числа e.

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

В конце концов, цель решения программистских задачек — не привнесение нового в науку, а поддержание мозгов в тонусе. Это и есть главная идея, лежащая в основе сайта "Зона кода".

Программа, решающая поставленную задачу, будет написана на языке C99. В ней будет использована библиотека big_int, в которой реализованы основные операции длинной арифметики. Но, для того, чтобы сконструировать алгоритм решения задачи, нам придётся прибегнуть к математическим выкладкам. Впрочем, для их понимания знакомства со стандартным двухгодичным курсом высшей математики, читаемым в технических вузах, хватит с лихвой.

Информация для проверки

Зайдя в Википедию и открыв статью, посвящённую числу e, я обнаружил, что в ней приведены 1000 знаков после запятой в десятичной записи этого числа (см. скриншот). Тогда я и решил ориентироваться на Википедию в плане количества находимых знаков, хотя до этого планировал остановиться на сотне цифр после запятой. Единственно, в чём я сомневался, так это в том, что мой компьютер справится с нахождением 1000 знаков за приемлемое время. Как оказалось, сомневался я не совсем уж напрасно.

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

Скриншот фрагмента статьи из Википедии, посвящённой числу e

В любом случае, сравнение нужно будет обязательно провести.

Предварительные замечания

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

Поясню, что библиотека big_int предназначена для работы с большими числами. Так я называю целые неотрицательные числа формально неограниченного размера. Но более распространённое название таких чисел — длинные целые. А операции, совершаемые с такими числами, принято называть длинной арифметикой.

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

Я буду предполагать, в дальнейшем, что читатель имеет хотя бы приблизительное представление о структуре big_int , описанной в библиотеке (именно в переменных типа big_int и хранятся большие числа), и о библиотечных функциях.

Что нам понадобится?

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

Ключевая формула, которую мы будем использовать, — это хорошо известное равенство числа e сумме числового ряда:

Этот ряд не просто сходится, но, что очень важно для нас, сходится достаточно быстро.

А ещё нам пригодится разложение функции e x по формуле Маклорена с остаточным членом в форме Лагранжа:

e x = ∑ n = 0 m x n n ! + e θ x x m + 1 m + 1 ! ,

где θ — это число, принадлежащее интервалу (0, 1). Сразу же перепишем это равенство для случая x = 1:

e = ∑ n = 0 m 1 n ! + e θ m + 1 ! .

Учитывая, что e θ < e < 3, получаем из последнего равенства следующую оценку:

e < ∑ n = 0 m 1 n ! + 3 m + 1 ! .

Оценка достаточно грубая, но нас она вполне устроит.

Представим в последнем неравенстве e в виде суммы ряда ∑ n = 0 ∞ 1 / n ! и вычтем из обеих частей неравенства сумму ∑ n = 0 m 1 / n ! :

∑ n = 0 ∞ 1 n ! - ∑ n = 0 m 1 n ! < 3 m + 1 ! .

В итоге, получаем:

∑ n = m + 1 ∞ 1 n ! < 3 m + 1 ! ,

Эта оценка будет нами использована в дальнейшем.

Построение алгоритма

Чтобы перевести задачу в "целочисленную плоскость", будем искать приближённое значение не e, а произведения 10 k e, где k — это некоторое натуральное число (о нём поговорим чуть позже). Зададимся целью получить оценку числа 10 k e вида

где p1 и p2 — это некоторые натуральные (k + 1)-разрядные числа, начинающиеся с цифры 2. Как только эта оценка будет получена, мы сравним эти два числа. Все идущие подряд цифры в записи числа p1, начиная со второй, совпадающие с соответствующими цифрами в записи p2, и будут, очевидно, являться первыми знаками после запятой в записи числа e. Нам нужно получить такую оценку, при которой количество таких цифр (обозначим его, например, s) будет не меньше 1000.

Ясно, что s зависит от k. Число k можно подбирать эмпирически. Однако, так как записи каждого из чисел p1 и p2 состоят из k + 1 разрядов, но, при этом, p1 ≠ p2, то s не может превышать k − 1. А поскольку мы должны получить s = 1000, то можно уже сейчас сказать, что k должно быть, по крайней мере, больше 1000.

Домножим теперь обе части равенства e = ∑ n = 0 ∞ 1 / n ! на 10 k . Имеем:

10 k e = ∑ n = 0 ∞ 10 k n ! .

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

p = ∑ n = 0 ∞ 10 k n ! .

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

Заметим, что, поскольку выражение 10 k / n! с возрастанием n убывает и стремится к 0 при n → ∞ , то, начиная с некоторого значения n, все числа вида 10 k / n! будут меньше единицы, т. е. их целые части будут нулевыми.

Другими словами, найдётся такой номер m, что будет выполнено:

10 k m ! ≥ 1 , но 10 k n ! < 1 при n > m , т.е. 10 k n ! = 0 при n > m .

Это означает, что, определив m таким образом, p теперь можно записать в виде:

p = ∑ n = 0 m 10 k n ! .

Найдём 10 k e − p:

10 k e - p = ∑ n = 0 ∞ 10 k n ! - ∑ n = 0 m 10 k n ! = ∑ n = 0 m 10 k n ! + ∑ n = m + 1 ∞ 10 k n ! .

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

Оценим разность 10 k e − p. Для этого сначала получим оценки для каждого из слагаемых, расположенных в правой части приведённого выше равенства. Начнём с первого из них. Заметим, что каждое из первых трёх слагаемых рассматриваемой нами суммы равно 0:

∑ n = 0 m 10 k n ! = 10 k 0 ! + 10 k 1 ! + 10 k 2 ! + ∑ n = 3 m 10 k n ! = 10 k + 10 k + 5 · 10 k - 1 + ∑ n = 3 m 10 k n ! = = ∑ n = 3 m 10 k n ! .

Заметим теперь, что

10 k 3 ! + 10 k 5 ! = 10 k 3 ! 1 + 1 4 · 5 = 10 k 6 · 21 20 = 7 · 10 k 40 = 175 · 10 k - 3 .

Но если сумма двух положительных дробных чисел целая, то сумма их дробных частей равна 1, т. е.

10 k 3 ! + 10 k 5 ! = 1 .

Возвращаемся к нашей сумме:

∑ n = 3 m 10 k n ! = 10 k 3 ! + 10 k 4 ! + 10 k 5 ! + ∑ n = 6 m 10 k n ! = 1 + 10 k 4 ! + ∑ n = 6 m 10 k n ! .

Оценим выражение, стоящее справа. Оно представляет собой сумму m − 3 слагаемых, причём, одно из них равно 1, а каждое из остальных строго меньше 1. Это означает, что вся сумма, которая, напомню, изначально была записана как ∑ n = 0 m 10 k / n ! , строго меньше m − 3:

∑ n = 0 m 10 k n ! < m - 3 .

Итак, оценка для первого слагаемого из суммы, в виде которой мы представили разность 10 k e − p, получена.

Переходим ко второму слагаемому. Вернёмся к полученной ранее оценке:

∑ n = m + 1 ∞ 1 n ! < 3 m + 1 ! .

Домножим обе части неравенства на 10 k :

∑ n = m + 1 ∞ 10 k n ! < 3 · 10 k m + 1 ! .

Учитывая, что, в силу выбора m, 10 k / (m + 1)! < 1, получаем из этого неравенства:

∑ n = m + 1 ∞ 10 k n ! < 3 .

Оценка для второго слагаемого получена! Сложим это неравенство с неравенством, представляющим собой оценку для первого слагаемого:

∑ n = 0 m 10 k n ! + ∑ n = m + 1 ∞ 10 k n ! < m + 3 - 3 = m .

Напомню, что сумма, стоящая слева, равна 10 k e − p. Заметим также, что эта сумма положительна. Таким образом,

0 < 10 k e − p < m, т. е. p < 10 k e < m + p.

Итак, искомая оценка для числа 10 k получена!

А теперь давайте сформулируем наш алгоритм кратко и по шагам.

Шаг 1. Выбираем число k > 1000.
Шаг 2. По заданному k находим числа p и m, такие, что

10 k m ! ≥ 1 , но 10 k m + 1 ! < 1 ,

В этом алгоритме ничего не сказано о том, какое число взять в качестве начального значение k и какими шагами увеличивать k в случае неудачных попыток. Об этом мы поговорим в следующем разделе.

Приближённое вычисление m и выбор начального значения k

Зададимся вопросом: сколько, хотя бы примерно, понадобится слагаемых для вычисления p? Как мы знаем, это число слагаемых превышает m на 1 (а всего слагаемых в ходе работы программы, включая последнее нулевое, будет вычислено m + 2). Давайте попробуем найти приближённое значение m.

Числа m и k, как известно, должны удовлетворять неравенствам

m! ≤ 10 k < (m + 1)!

Другими словами, m должно быть таково, чтобы число 10 k попадало в промежуток [m!, (m + 1)!) . Чтобы найти примерное значение m, можно приближённо приравнять одну из границ промежутка (например, левую), числу 10 k :

m! ≈ 10 k .

Давайте теперь приблизим m! с помощью формулы Стирлинга, достаточно точной для больших значений m (превышающих 100):

m ! ≈ 2 π m m e m .

Подставим m! из последнего приближённого равенства в предпоследнее:

2 π m m e m ≈ 10 k .

Возьмём логарифм по основанию 10 от обеих частей этого приближённого равенства:

m lg m - lg e + lg m + lg 2 π 2 ≈ k .

Наша задача — решить это уравнение относительно m. Будем решать его графически. Для этого построим график функции

y = x lg x - lg e + lg x + lg 2 π 2 .

Для построения графика я использовал Mathcad. График приведён ниже.

График функции y = x (lg x − lg e) + (lg x + lg (2 π)) / 2

Как мы видим, на отрезке [400, 500] график неотличим от отрезка прямой. Кстати, второе слагаемое в функции (дробь с двойкой в знаменателе) даёт очень несущественный вклад на этом отрезке. График функции без второго слагаемого крайне близок к графику исходной функции.

Вернёмся к нахождению приближённого значения m. Для начала нам нужно определиться со значением k. Мы знаем лишь, что k немного превышает 1000. Положим для простоты: k = 1000. Тогда для нахождения m проведём прямую линию, параллельную оси абсцисс, через точку с координатами (0, 1000) до её пересечения с графиком функции, после чего опустим перпендикуляр на ось абсцисс (на рисунке эти две линии обозначены синим пунктиром). Основание перпендикуляра попадает в точку с координатами (450, 0).

Таким образом, функция принимает значение, равное 1000, в точке x ≈ 450, откуда m ≈ 450. Приближённое значение m получено!

Количество слагаемых, сумма которых равна p, тоже можно считать примерно равным 450. Что ж, это вполне приемлемое для нас количество. Обычный компьютер должен справиться с вычислениями, потратив на это не слишком много времени.

Поскольку число m — трёхзначное, то 3 последние цифры у чисел p и p + m совпадать точно не будут, поэтому значение s не будет превышать k − 3. Таким образом, k не может быть меньше 1003. Именно 1003 и возьмём в качестве начального значения k. Далее, при необходимости, будем увеличивать значения k с шагом 1.

А теперь переходим к программированию!

Создание программы

Создадим файл calc_e.c и подключим к нему заголовочный файл библиотеки big_int:

Программу создадим в двух версиях. Первая версия — тестовая. Она выведет на консоль числа k, m, s, вычислив, предварительно, второе и третье из них. Увеличивая k, будем добиваться того, чтобы s достигло 1000. Как только k, гарантирующее это достижение, будет получено, запустим вторую версию — рабочую. Она, используя найденное значение k, вычислит нужные нам знаки в десятичной записи числа e и напечатает их.

Управлять версиями программы будем через макрос TEST :

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

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

Помимо приведённых строк, файл calc_e.c будет содержать только функцию main() , которая и будет "отвечать" за все вычисления, описанные в этой статье. Ниже приведён код этой функции.

Перед тем, как прокомментировать код функции main() , опишу смысл всех задействованных в ней переменных типа big_int , идентифицируя их с помощью указателей.

tenpk (сокращение от ten power k) указывает на переменную, содержащую число 10 k .
p указывает на переменную, предназначенную для хранения частичных сумм ряда, сходящегося к p (напомню, что только первые m + 1 членов ряда отличны от нуля). Окончательное значение переменной совпадает с p.
m указывает на переменную, предназначенную для хранения уменьшенных на единицу номеров слагаемых, участвующих в вычислении частичных сумм ряда, сходящегося к p.
fact указывает на переменную, предназначенную для хранения факториалов чисел, хранимых в переменной, адресуемой m .
one указывает на переменную, содержащую число 1.
temp указывает на переменную, содержащую слагаемое, участвующее в вычислении частичных сумм ряда, сходящегося к p.
ppk (сокращение от p plus k) указывает на переменную, содержащую сумму p + k.

Работа функции main() начинается с создания символьного массива str (стр. 3), в который помещается текстовое представление (в виде C-строки) числа 10 k (см. стр. 4-6). Далее на основе созданной строки формируем переменную, содержащую 10 k (стр. 7). Значения её полей сохраняем в отдельных переменных (стр. 8, 9). Нам они пригодятся в дальнейшем.

Создаём переменную для хранения частичных сумм ряда, сходящегося к p. Сумму первых двух членов ряда не будем вычислять с помощью цикла, а сразу же поместим в эту переменную (стр. 10). Как видите, её значение получено сложением содержимого переменной, адресуемой указателем tenpk , с ним самим.

Переменным, адресуемым указателями fact и m , при создании присваиваем единичные начальные значения (стр. 11, 12). Создаём и инициализируем переменную, содержащую единицу (стр. 13) и создаём указатель temp (стр. 14).

Далее следует основная часть функции main() — цикл while , в ходе которого и вычисляется сумма, равная p. Сами слагаемые вычисляются, начиная с 3-го, а сумма "накапливается" в переменной, адресуемой p .

Давайте разберём условие продолжение цикла while (см. стр. 15, 16).

Создаётся копия переменной, содержащей число 10 k , с помощью функции create() . Значение переменной, адресуемой m , увеличивается на 1 посредством функции add() . Полученный результат перемножается с использованием функции multiply() с содержимым переменной, предназначенной для хранения факториалов , и новое значение факториала сохраняется в этой же переменной. Наконец, переменная, содержащая 10 k , делится нацело на новый факториал с помощью функции divide() , результат сохраняется в этой же переменной, а её адрес помещается в temp . Значение старшего байта результирующей переменной и является условием продолжения цикла.

Резюмируем предыдущий абзац: в результате всех этих действий создаётся переменная, содержащая текущее слагаемое. Кроме того, обновляются значения переменных, адресуемых указателями m и fact . Старший байт созданной переменной будет отличен от нуля тогда и только тогда, когда само число, хранящееся в ней, не равно 0. В этом случае цикл продолжится.

В теле цикла добавляем вычисленное слагаемое к частичной сумме, хранящейся в переменной, адресуемой p (результат сложения сохраняется в этой же переменной) (стр. 18). Удаляем найденное слагаемое за ненадобностью (стр. 19).

Как только очередное найденное слагаемое окажется нулевым, цикл завершится. В переменной, адресуемой p , после завершения цикла будет находиться искомый результат — число p. А в переменной, адресуемой m , — число m + 1.

Удаляем все уже ненужные нам переменные (стр. 21-23) и создаём строку, содержащую текстовое представление p (стр. 24).

Далее идёт код, выполняющийся только в тестовом режиме. Вычисляем сумму p + m и помещаем её в новую переменную (стр. 26), после чего получаем её текстовое представление (стр. 27), а саму переменную удаляем (стр. 28).

Заметим, что в строке 26, помимо прочего, был выполнен декремент переменной, адресуемой указателем m , в результате чего её значение сравнялось с числом m.

Далее, анализируя текстовые строки, содержащих десятичные записи чисел p и p + m, находим число s и помещаем его в переменную s (стр. 29-32). Удаляем текстовое представление p + m (стр. 33). Наконец, выводим на консоль числа k, m и s (стр. 34).

Обратите внимание на 3-й аргумент функции printf() (стр. 34). Мне не хотелось вызывать функцию, возвращающую текстовое представление переменной, адресуемой m , поэтому я решил прибегнуть к трюку. Мы знаем, что число m близко к 450, поэтому для его хранения будет использоваться двухбайтовый массив. Этот массив можно интерпретировать как ячейку памяти, в которой хранится беззнаковое короткое целое число. Чтобы получить его, достаточно привести адрес первого элемента массива к типу unsigned short* и разыменовать полученный адрес, что и было проделано.

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

Мы выводим на консоль текст "e = 2,", а затем печатаем искомые 1000 десятичных знаков после запятой, завершающиеся многоточием (см. стр. 36-50). Я решил прибегнуть к тому же способу размещения знаков, который использован в Википедии: 50 знаков в каждой строке, разбитых пробелами на 5 групп по 10 знаков в каждой. Так проще будет сравнивать наш результат с приведённым в Википедии.

Далее идёт завершающий фрагмент функции main() , выполняющийся в обоих режимах. Удаляем оставшиеся неудалёнными переменные типа big_int (стр. 52-54) и строку, содержащую текстовое представление p (стр. 55).

Выполнение программы

Когда я запустил программу впервые, она, как мне показалось, очень долго не хотела завершаться. Я уже начал переживать: а вдруг, подумал я, для вычислений потребуется время, равное паре продолжительностей жизни вселенной? Но нет, 30-ти миллиардов лет не понадобилось. Мне показалось, что программа работала минут 5, но потом оказалось, что гораздо меньше. Причина такого эффекта крылась, видимо, в моём нетерпении, повлиявшем на восприятие времени.

Итак, напомню, что сначала мы используем тестовый режим и значение k, равное 1003. Вот что выводится на консоль:

Ну что ж, уже неплохой результат! Не хватает единицы до полного счастья. Заметим также, что приближённое значение m было найдено достаточно точно. Кстати, мы находили его для k = 1000 и получили, напомню, 450. А если бы искали для k = 1003, то 451 и получили бы.

Увеличиваем k на 1, изменяя соответствующую строку кода:

Компилируем и запускаем программу снова. Вот результат:

Изменилось только k, и больше ничего. Опять увеличиваем k:

Вот что получаем:

Раньше был недобор, а теперь — перебор. Но нас он вполне устраивает! Итак, останавливаемся на k = 1005.

Переключаемся с тестового режима на рабочий, превращая в комментарий строку, определяющую макрос TEST:

Компилируем и запускаем. Получаем следующий консольный вывод:

Заключение

Десятичные знаки после запятой в записи числа e, полученные нами, полностью совпадают с приведёнными в Википедии, поэтому в их достоверности можно не сомневаться.

Способ решения задачи, выбранный нами, помог нам убедиться в том, что с помощью целочисленной арифметики иногда можно решать задачи, в которых фигурируют вещественные числа. Кстати, в нашей программе использовались все четыре арифметических операции, реализованные в big_int (правда, возможность получения остатка от целочисленного деления оказалась невостребованной). Можно сказать, что библиотека big_int с честью выдержала ещё одну "проверку на прочность".

Я создавал и запускал программу на своём стареньком, видавшем виде компьютере, управляемом операционной системой Windows XP. Установленный на нём процессор AMD Athlon 64 3400+, работающий на частоте 2,40 ГГц, считался очень неплохим лет 10 назад. Объём оперативной памяти компьютера, равный 1 ГБ, тоже, по тем временам, был вполне приемлемым для комфортной офисной работы.

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

Впрочем, "обкатывать" изменения можно, устанавливая значительно более скромные начальные данные, чем требуется для решения задачи. К слову сказать, 100 знаков после запятой у числа e программа находит менее, чем за секунду.

Ну, а естественным продолжением поднятой темы было бы нахождение 1000 знаков после запятой в десятичной записи числа π. Уверен, что и до решения этой задачи руки обязательно дойдут.

📎📎📎📎📎📎📎📎📎📎