Формула Бэйли-Боруэйна-Плаффа

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

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

Сама формула выглядит следующим образом:

А так выглядит реализация, которую мы оставляем без подробных разъяснений (при этом требуется лишь ряд параметров, в частности, позиция цифры в числе Пи и сколько цифр с этой позиции вычислить):

import std.conv;
import std.math;
import std.stdio;

auto ihex(double x, int nhx)
{
	enum string HEX = "0123456789ABCDEF";
	char[16] result;
	double y = fabs(x);

	for (int i = 0; i < nhx; i++)
	{
		y = 16.0 * (y - floor(y));
		result[i] = HEX[to!int(y)];
	}

	return result;
}

double series(int m, int id)
{
	double eps = 1e-17, s = 0.0;
	double ak, p, t;
	
	for (int k = 0; k < id; k++)
	{
		ak = 8 * k + m;
		p = id - k;
		t = expm(p, ak);
		s +=t / ak;
		s -= to!int(s);
	}

	for (int k = id; k <= id + 100; k++)
	{
		ak = 8 * k + m;
		t = (16.0 ^^ to!double(id - k)) / ak;
		
		if (t < eps)
		{
			break;
		}

		s += t;
		s -= to!int(s);
	}
	return s;
}


double expm(double p, double ak)
{
	double p1, pt, r;
	enum int NTP = 25;
	static double[NTP] tp;
	static int tp1 = 0;
	int i;

	if (tp1 == 0)
	{
		tp1 = 1;
		tp[0] = 1;

		for (i = 1; i < NTP; i++)
		{
			tp[i] = 2.0 * tp[i - 1];
		}
	} 

	if (ak == 1.0)
	{
		return 0.0;
	}

	for (i = 0; i < NTP; i++)
	{
		if (tp[i] > p)
		{
			break;
		}
	}

	pt = tp[i - 1];
	p1 = p;
	r = 1;

	for (int j = 1; j <= i; j++)
	{
		if (p1 >= pt)
		{
			r = 16.0 * r;
			r = r - to!int(r / ak) * ak;
			p1 = p1 - pt;
		}

		pt = 0.5 * pt;

		if (pt >= 1.0)
		{
			r *= r;
			r -= to!int(r / ak) * ak;
		}
	}

	return r;
}

void main()
{
        // номер позиции
	int ID = 5;
        // количество цифр
	int NHX = 16;
	char[] chx;

	auto s1 = series(1, ID);
	auto s2 = series(4, ID);
	auto s3 = series(5, ID);
	auto s4 = series(6, ID);
	
	double PID = 4.0 * s1 - 2.0 * s2 - s3 - s4;
	PID = PID - to!int(PID) + 1.0;
	
	writefln(" position = %d\n fraction = %.15f\n hex digits = %10.10s",
		ID,
		PID,
		ihex(PID, NHX)
		);
}

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