abiao 发表于 2010-8-12 19:26:27

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

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

gxqcn 发表于 2010-8-13 09:11:52

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

mathe 发表于 2010-8-13 10:36:05

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

liangbch 发表于 2010-8-13 11:44:54

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

abiao 发表于 2010-8-13 19:38:16

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

abiao 发表于 2010-8-14 17:46:41

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

〇〇 发表于 2010-8-15 16:07:39

源代码呢

abiao 发表于 2010-8-15 16:32:31

源代码比较多,贴出乘法的吧
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 = Num_Array;
      if(cutOff == 0 && Num_Array != 0)
      {
            cutOff = Precision - i;
      }
      Num_Array = 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;

      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;
                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 // low:eax,high:edx
                        add eax,
                        adc edx,0
                        add eax,ebx
                        adc edx,0
                        mov ,eax
                        mov ebx,edx
                  loop lp
                  mov carry,ebx
                }
            }

            //继续进位
            if (shift>0)
            {
                unsigned int* p0 = &Num_Array;
                _asm
                {
                  mov ecx,shift
                  mov edi,p0
                  mov edx,carry
                  lp2:
                        mov eax,
                        add eax,edx
                        mov ,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 = Num_Array;
                }
                Num_Array = carry;

            }
      }
    }

    (*pExponent) += temp_Exponent - 1;
}

abiao 发表于 2010-8-16 19:48:57

我用了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 发表于 2010-8-17 01:07:26

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

就楼主的算法而言,每次迭代精度加倍,应该是比较高效的算法了。麻烦的在于计算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)时收敛的很慢。
页: [1] 2 3 4 5
查看完整版本: 我想的圆周率计算法,1万位十进制精度计算时间在数秒至数十秒左右