Миллион цифр числа "пи"
Миллион цифр числа "пи"
Хотите готовиться со мной к ЕГЭ?
Пишите: ydkras@mail.ru
Немного обо мне.
(Данная тема имеет довольно слабое отношение к ЕГЭ по информатике, но оказалась красивой.)
Полвека назад (страшно сказать!) я участвовал в институтской олимпиаде по программированию. Одной из задач было вычислить, какая цифра стоит в числе "пи" то ли на 50-м, то ли на 100-м месте после запятой.
Я использовал ряд для арктангенса. Понятно, 4*arctg(1)=Pi. Разумеется, точности представления вещественных чисел в ЭВМ БЭСМ-4М было совершенно недостаточно, чтобы понять, что там творится в 50-м знаке после запятой. Поэтому пришлось, само собой, представлять числа в виде массивов цифр и реализовать процедуры для арифметических действий с такими числами - на языке Алгол-60 (многие ли читатели слышали про такой, интересно?).
Наконец, программа (на перфокартах, понятное дело - впрочем, мне приходилось иметь дело и с перфолентами) была отлажена и стала выдавать цифру за цифрой (очень неспешно). В итоге я занял на той олимпиаде 3-е место и получил денежный приз в размере примерно половины месячной стипендии.
Полвека пролетело как-то незаметно, и несколько дней назад я увидел, что некий студент ищет помощи: в институте ему дали задание вычислить 1000000 (прописью - один миллион) цифр числа "пи". С одной стороны - если человек просит помощи, ему надо помочь. С другой - я вспомнил, как я решал тогда задачу и сколько времени я на нее потратил, и желание помогать почти угасло.
Но тут в готову мне пришла простая мысль: если студентам дают такое задание, то... оно решаемо? Но как? Чтобы получить миллион верных цифр, надо суммировать ну очень много членов ряда. Ряд нужен какой-то другой, быстро сходящийся. Опять же надо писать реализацию арифметических операций над числами-массивами. Ну хорошо, в Питоне можно работать с целыми числами из сколь угодно большого количества цифр, но мне-то нужны дробные.
На помошь пришел, естественно, гугл. Я нашел очень интересную статью про число "пи" и методы его вычисления: "Число пи". Из неё я узнал, что в то время, когда я решал задачу на олимпиаде, никто в мире ещё не получил миллион цифр числа "пи". Однако с тех пор математики существенно продвинули эту тему (да и компьютеры стали несколько порезвее, чего уж там). В частности, был открыт алгоритм Брента-Саламина, который с первого взгляда мне понравился своей компактностью.
Алгоритм Брента-Саламина
Начальные установки.
Положим
a0=1
b0=1/sqrt(2)
t0=1/4
p0=1
Цикл алгоритма.
an+1=(an+bn)/2
bn+1=sqrt(an*bn)
tn+1=tn-pn*(an-an+1)**2
pn+1=2*pn
Цикл продолжается, пока числа an и bn не станут достаточно близки друг другу. Оценка числа "пи" дается формулой
pi=(an+bn)**2/(4*tn)
Утверждалось, что 25 итераций цикла достаточно, чтобы получить 45 миллионов (!) верных цифр. Мне нужно в 45 раз меньше. Если на каждом шаге число верных цифр удваивается (как утверждалось в упомянутой выше статье), то мне нужно всего-то примерно 20 итераций. Но... арфиметики для этого алгоритма мало, надо ещё и квадратный корень уметь извлекать! Беда.
Ну хорошо, можно написать процедуру итеративного вычисления корня. Но нельзя ли как-то попроще? А что может предложить питон? И снова я обратился к гуглу.
Оказалось, что в питоне можно работать не только с вещественноми числами (float), но и с десятичными дробями (decimal). Более того, для величин decimal реализованы квадратный корень, экспонента и натуральный логарифм. Это плюс. Минус же состоял в том, что эти десятичные числа состоят из... 28 десятичных цифр - а мне нужен миллион! Но оказалось, что длину чисел по умолчанию можно менять, и питон не имеет никаких возражений против чисел decimal из миллиона и более цифр!
Более подробно про числа decimal и работу с ними можно прочитать здесь: https://stepik.org/lesson/360941/step/1?unit=345464
В итоге всё вылилось в вот такую программу из менее чем двух десятков строк:
from decimal import *
getcontext().prec = 1100000
a0=Decimal("1")
b0=1/Decimal("2").sqrt()
t0=Decimal("0.25")
p0=Decimal("1")
pi0=3
while True:
a1=(a0+b0)/2
b1=(a0*b0).sqrt()
t1=t0-p0*(a0-a1)**2
p1=2*p0
pi1=(a1+b1)**2/(4*t1)
if pi0==pi1:
break
a0,b0,t0,p0,pi0=a1,b1,t1,p1,pi1
pi=str(pi1)[:1000002]
# Сохраните строку pi в файле - или напечатайте, если хотите :)
Итак,час моего времени и почти пять минут времени моего ноутбука - и миллион цифр числа "пи" получен. Как и предполагалось, для этого потребовалось 20 итераций, каждая итерация требовала примерно 13 секунд.
Несколько заключительных замечаний.
1. Числа decimal - простой и удобный способ вычислений практически с любой желаемой точностью - буквально в миллионы цифр после запятой, как мы имели возможность убедиться. Однако по скорости этот способ гораздо хуже, чем использование типа float.
2. В питоне есть также и рациональные дроби (т.е. числа вида m/n, где m - целое число, а n - натуральное). Для работы с ними надо импортировать модуль fraction. Кратко познакомиться с этой темой можно здесь: "Модуль fractions".А здесь - интересный обзор числовых типов в питоне: https://habr.com/ru/companies/wunderfund/articles/647331/
3. Комплексные числа в питоне тоже есть. Очень полезно для электротехнических расчетов, например. Подробнее смотри, например. тут: https://tirinox.ru/complex-python-mandelbrot/
4. Приведенная выше программа - любопытный тест для оценки производительности компьютера. Как уже было сказано, на моем ноутбуке (уже не первой молодости) она выполняется менее чем за 5 минут. А на стареньком Macbook 2006 года выпуска она трудилась почти 13 минут.
Хочется закончить этот рассказ цитатой из статьи "Трансцендентное число "Пи"" глубоко уважаемого мной Мартина Гарднера (взято отсюда: https://coollib.net/b/147528/read#t42, рекомендую к прочтению):
"Вычисление нескольких тысяч знаков "пи" в настоящее время стало популярным средством проверки новых компьютеров и обучения молодых программистов. «Загадочное и чудесное "пи",— пишет в своей книге «Что мы знаем о больших числах» Филипп Дж. Дэвис, — стало чем-то вроде покашливания, которым компьютеры прочищают горло»."
(c) Ю.Д.Красильников, 2023 г.
Комментарии
Отправить комментарий