我想的圆周率计算法，１万位十进制精度计算时间在数秒至数十秒左右

abiao

我用自制的长浮点数计算ｐｉ，１万位十进制精度计算时间在数秒至数十秒左右，视机器性能而定，
基本算法是
x=3或者3.1或者3.14等等;
循环(x+=sin(x));
每计算一次sin可以获得一倍的精度，sin是用级数叠加的，每叠加一次需要计算4次乘法，一次除法，一次加法，获得２位左右的10进制精度
加乘除法是用最常规的方法计算的，没有使用傅立叶变换，只是在汇编级进行了优化，利用了x86的指令集的特性
基本的数据结构是使用一个无符号整数数组，每个都是一个计算单位,(2^32=4G)进制，表示时再转回10进制，一个符号标志，一个指数标志
性能在使用常规方法计算中应该算快的了吧
这个计算方法是在我算sin(3.14)=0.00159xxxx时想到的

gxqcn

我估计这样的算法流程是比较慢的。
因为：
1、计算pi的数学公式不够高效；
2、未使用高级大数乘法算法（FFT之类）；
3、先2进制，再转10进制，这个转换是非常耗时的，如果处理不好的话。

mathe

最主要是sin函数本身效率不会很高。
至于说高级大数乘法运算，是具体如何做乘法的问题，倒是同楼主这里的算法没有关系。
二进制到十进制转化倒还好，这个转化过程如果实现得好，相对于计算过程花费的时间是很少的。我记得最近计算派的记录的那个链接里面就是采用2进制计算的，然后最后转化为10进制的

liangbch

计算sin(x)本身就很麻烦，这样会很慢。常规算法中，Martin公式应该比你的算法更快，而且更简单。

abiao

算pi的主要目的是测试三角函数的正确性

abiao

另外Martin公式太复杂，我的这个比较简洁，虽然说sin要计算多次，不过前几次的时间可以忽略不计的

〇〇

源代码呢

abiao

源代码比较多，贴出乘法的吧
void mul(int Precision, unsigned int* Num_Array, int* pSign, int* pExponent, unsigned int* num_Num_Array, int num_Sign, int num_Exponent, unsigned int* buf_Num_Array)
{
    int i;
   
    int cutOff = 0;

    for (i = Precision - 1; i >=0; i--)
    {
        buf_Num_Array[i] = Num_Array[i];
        if(cutOff == 0 && Num_Array[i] != 0)
        {
            cutOff = Precision - i;
        }
        Num_Array[i] = 0;
    }
    cutOff--;

    int temp_Exponent = (*pExponent) + num_Exponent;
    *pExponent = 0;

    *pSign = (*pSign) * num_Sign;

    for (int idxMain = 0; idxMain < Precision; idxMain++)
    {
        int n = num_Num_Array[idxMain];

        if (n!=0)
        {
            int shift = idxMain + (*pExponent);//(*pExponent)修正量
            int cnt = Precision - (cutOff > shift ? cutOff : shift);

            int carry = 0;
            if (cnt>0)
            {
                unsigned int* p0 = &Num_Array[shift - 1];
                unsigned int* p1 = &buf_Num_Array[-1];
                _asm
                {
                    mov ecx,cnt
                    mov edi,p0
                    mov esi,p1
                    mov ebx,0
                    lp:
                        mov eax,n
                        cdq
                        mul DWORD PTR [esi+ecx*4] // low:eax,high:edx
                        add eax,[edi+ecx*4]
                        adc edx,0
                        add eax,ebx
                        adc edx,0
                        mov [edi+ecx*4],eax
                        mov ebx,edx
                    loop lp
                    mov carry,ebx
                }
            }

            //继续进位
            if (shift>0)
            {
                unsigned int* p0 = &Num_Array[shift - 1];
                _asm
                {
                    mov ecx,shift
                    mov edi,p0
                    mov edx,carry
                    lp2:
                        mov eax,[edi]
                        add eax,edx
                        mov [edi],eax
                        mov edx,0
                        jnc ext
                        mov edx,1
                        sub edi,4
                    loop lp2
                    ext:
                    mov carry,edx
                }
            }

            if (carry!=0)
            {
                (*pExponent)++;
                for (i = Precision - 1; i > 0; i--)
                {
                    Num_Array[i] = Num_Array[i - 1];
                }
                Num_Array[0] = carry;

            }
        }
    }

    (*pExponent) += temp_Exponent - 1;
}

abiao

我用了Martin公式(pi=16*arctan(1/5)-4*arctan(1/239))
发现比我的算法还慢

因为sin的收敛的速度比atan快，叠加的次数少
sin x = x-x^3/3!+x^5/5!-...(-1)k-1*x2k-1/(2k-1)!+... (-∞<x<∞)
arctan x = x - x^3/3 + x^5/5 - ... (x≤1)
sin式的被除数部分是阶乘，arctan是累加（？）

实际计算时，sin的x是3.14xxx，arctan是1/5,1/239
即使如此还是sin比较快

liangbch

我试了下，楼主的算法确实可行。

就楼主的算法而言，每次迭代精度加倍，应该是比较高效的算法了。麻烦的在于计算sin(x)
  计算sin(x)是不需要计算出全部精度的，比如当x=3时，计算sin(x)只需要精确2位小数，当x=3.14时，计算sin(x)只需要精确4位小数，当x越接近于pi,需要的sin(x)的精度越高。
  但是，这里需要解决一个核心问题，如何计算sin(x)，当x接近于pi时，用级数展开式计算sin(x)时收敛的很慢。