Я реализовал метод Bostan–Gaudry–Schost для хвостовых сумм. По числу операций это Õ(n^{3/2}) вместо O(n²). Но по времени он начинает выигрывать только при n порядка 10⁷, а на всех n, которые я реально прогонял, он проигрывает прямому счёту. Скрипт — pi_decimal_digits.py: если хвост длиннее порога THR (по умолчанию 150000), его считает метод BGS. Порог можно задать четвёртым аргументом командной строки.
Что сделано математически
Остаётся только неполная сумма на хвосте подблока, то есть префикс гипергеометрического члена по модулю a. Матрицы 2×2 по шагам перемножаются схемой BGS: произведение блоков длины v≈√(L/3) считается удвоением со сдвигом точек выборки по Лагранжу (свёртки по модулю a).
Что измерено
Корректность. Хвосты метода совпали с прямым счётом в 25 случайных тестах. Целиком программа с принудительным BGS (порог 200 и 300 вместо 150000) выдаёт верные цифры при n=16000 и n=8000, но заметно медленнее прямого счёта: 18 с против 1,6 с при n=16000.
Рост по длине хвоста L. Время одного хвоста растёт как √L: 0,029 с при L=10⁶ и 0,263 с при L=6,4·10⁷. Отношение t/√L остаётся ≈0,03 мс.