C# BBP公式计算π

BBP公式可以直接计算出小数点后d位的数,不依赖第d位之前的数字,计算结果为16进制.
BBP公式 :

1. 核心公式
   π = ∑ k = 0 ∞ { 1 1 6 k ( 4 8 k − 1 − 2 8 k + 4 − 1 8 k + 5 − 1 8 k + 6 ) } (1) \pi=\sum_{k=0}^{\infty}\{\frac{1}{16^{k}}(\frac{4}{8k-1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6})\}\tag{1} π=k=0∑∞​{16k1​(8k−14​−8k+42​−8k+51​−8k+61​)}(1)
2. 推导公式
   { 1 6 d π } = { 4 { 1 6 d S 1 } − 2 { 1 6 d S 4 } − { 1 6 d S 5 } − { 1 6 d S 6 } } (2) \{16^{d}\pi\}=\{4\{16^{d}S_{1}\}-2\{16^{d}S_{4}\}-\{16^{d}S_{5}\}-\{16^{d}S_{6}\}\}\tag{2} {16dπ}={4{16dS1​}−2{16dS4​}−{16dS5​}−{16dS6​}}(2)
   其中
   $$S_j=\sum _{k=0}^{\infty }\frac{1}{16^k(8k+j)}\text{\hspace{1em}(3)}$$
3. 推导整理出用于计算的公式
   { 1 6 d S j } = { { ∑ k = 0 d 1 6 d − k 8 k + j } + ∑ k = d + 1 ∞ 1 6 d − k 8 k + j } = { { ∑ k = 0 d 1 6 d − k   m o d   8 k + j 8 k + j } + ∑ k = d + 1 ∞ 1 6 d − k 8 k + j } (4) \begin{aligned} \{16^{d}S_{j}\}&=\{\{\sum_{k=0}^{d}\frac{16^{d-k}}{8k+j}\}+\sum_{k=d+1}^{\infty}\frac{16^{d-k}}{8k+j}\}\\ &=\{\{\sum_{k=0}^{d}\frac{16^{d-k}\bmod8k+j}{8k+j}\}+\sum_{k=d+1}^{\infty}\frac{16^{d-k}}{8k+j}\} \end{aligned}\tag{4} {16dSj​}​={{k=0∑d​8k+j16d−k​}+k=d+1∑∞​8k+j16d−k​}={{k=0∑d​8k+j16d−kmod8k+j​}+k=d+1∑∞​8k+j16d−k​}​(4)
   公式相关的说明戳这里
   按公式2码出来的. 上码

using System;
using System.Diagnostics;

namespace BBP
{
    public class BBP
    {
        /// <summary>
        /// 二元展开快速求解 POW. 减少乘法计算次数
        /// </summary>
        /// <param name="a">底数</param>
        /// <param name="b">指数</param>
        /// <returns>计算结果</returns>
        static long QuickPow(long a, long b)
        {
            long rlt = 1;
            while (b != 0)
            {
                if ((b & 1) != 0)
                {
                    rlt *= a;
                }
                a *= a;
                b /= 2;
            }
            return rlt;
        }

        /// <summary>
        /// 二元展开快速求解 POW & MOD
        /// </summary>
        /// <param name="a">底数</param>
        /// <param name="b">指数</param>
        /// <param name="mod">除数</param>
        /// <returns>计算结果</returns>
        static long QuickPowMod(long a, long b, long mod)
        {
            long rlt = 1;
            while (b != 0)
            {
                if ((b & 1) != 0)
                {
                    rlt = rlt * a % mod;
                }
                a = a * a % mod;
                b /= 2;
            }
            return rlt;
        }

        /// <summary>
        /// 计算 16^d * Sj
        /// </summary>
        /// <param name="d">第几位</param>
        /// <param name="j">1/4/5/6</param>
        /// <returns></returns>
        static double Calculate16dSj(int d, int j)
        {
            double part1 = 0;
            double part2 = 0;

            // 精度按 d ~ d + 10 
            int accuracy = d + 10;

            for (int k = 0; k <= d; k++)
            {
                part1 += QuickPowMod(16, d - k, 8 * k + j) * 1.0 / (8 * k + j);
            }

            for (int k = d + 1; k < accuracy; k++)
            {
                part2 += QuickPow(16, d - k) / (8 * k + j);
            }

            return part1 + part2;
        }

        /// <summary>
        /// 计算 16^d * PI
        /// </summary>
        /// <param name="d">第几位</param>
        /// <returns></returns>
        static double Calculate16dPI(int d)
        {
            double rlt = 4 * Calculate16dSj(d, 1) - 2 * Calculate16dSj(d, 4) - Calculate16dSj(d, 5) - Calculate16dSj(d, 6);

            // 只需要小数部分,所以减去整数部分
            rlt -= (int)rlt;

            // 因为在计算 16^d * Sj 的第一部分时用到了取模, 所以可能计算结果是负数, 需要加 1
            if (rlt < 0) rlt += 1;
            return rlt;
        }

        /// <summary>
        /// 按顺序逐位计算
        /// </summary>
        /// <param name="end">到第几位结束</param>
        /// <param name="start">从第几位开始,默认第一位(从 0 开始)</param>
        static void CalculateInOrder(int end, int start = 0)
        {
            Console.WriteLine($"开始计算 π 第 {start + 1} ~ {end} 位");
            for (int d = start; d < end; d++)
            {
                Console.WriteLine($"小数点后第 {d + 1} 位 -> {(int)(Calculate16dPI(d) * 16):x}");
            }
        }

        /// <summary>
        /// 只计算某一位, 第一位是从 0 开始
        /// </summary>
        /// <param name="index">第几位, 从 0 开始. 例如:想算第 3 位, index = 2</param>
        static void CalculateSpecifyIndex(int index)
        {
            Console.WriteLine($"开始计算 π 第 {index + 1} 位");
            Console.WriteLine($"小数点后第 {index + 1} 位 -> {(int)(Calculate16dPI(index) * 16):x}");
        }

        static void Main(string[] args)
        {
            Console.WriteLine("计算 π (十六进制)\n准备就绪...");
            Stopwatch stopwatch = new Stopwatch();
            stopwatch.Start();

            // 先算算前20位
            CalculateInOrder(20);
            // 再算算小数点后第1位
            CalculateSpecifyIndex(0);

            stopwatch.Stop();
            Console.WriteLine($"共用时 {stopwatch.ElapsedMilliseconds} ms");
        }
    }
}

EOF

参考文档:
圆周率pi的BBP计算公式之详细证明
BBP算法计算圆周率（BBP Formula HDU - 6217）