largefloat.cpp

一个很叨 的算法作业,很强大的源代码 很强大的源代码很强大的源代码

///////////////////////////////////////////////////////////////////////////////////
// LargeFloat.cpp
//

#include <iostream>
#include "TLargeFloat.h"

//#define _MY_TIME_CLOCK

#ifdef _MY_TIME_CLOCK
    #include <windows.h>
    #define clock_t __int64
    clock_t CLOCKS_PER_SEC=0;
    inline clock_t clock() 
    {
        __int64 result;
        if (CLOCKS_PER_SEC==0)
        {
            QueryPerformanceFrequency((LARGE_INTEGER *)&result);
            CLOCKS_PER_SEC=(clock_t)result;
        }
        QueryPerformanceCounter((LARGE_INTEGER *)&result);
        return (clock_t)result;
    }
#else
    #include <time.h>
#endif


//计算圆周率PI
TLargeFloat GetBorweinPI(unsigned long uiDigitsLength);
TLargeFloat GetAGMPI(unsigned long uiDigitsLength);
TLargeFloat GetChudnovskyPI(unsigned long uiDigitsLength);

void Debug_toCout(const std::string& strx,const TLargeFloat& x)//调试输出
            { std::cout<<strx<<std::endl<<x<<std::endl; }

int main()
{
    try
    {
        LargeFloat_UnitTest();

        TLargeFloat x(0.0,TLargeFloat::TDigits(100));
        x=1;
        Debug_toCout("1/7:=",x/7);

        x.StrToLargeFloat("12345678987654321");
        Debug_toCout("12345678987654321*12345678987654321 := ",x*x);

        x=2;
        Debug_toCout("2^0.5 := ",sqrt(x));

        //test PI
        std::cout<<std::endl<<"请输入需要PI的有效位数:";
        long pi_length=0; std::cin>>pi_length;    std::cout<<std::endl;
        clock_t t0=clock();
        x=GetAGMPI(pi_length); //////////////////////////////////
        //x=GetBorweinPI(pi_length) 
        t0=clock()-t0;

        Debug_toCout("PI := ",x);
        std::cout<<"PI计算时间:"<<(t0*1.0/CLOCKS_PER_SEC)<<" 秒";

    }
    catch (const TLargeFloatException& lfe)
    {
        std::cout<<">> LargeFloat运算时发生异常: "<<lfe.what()<<std::endl;
    }
    catch (const std::exception& e)
    {
        std::cout<<">> 异常: "<<e.what()<<std::endl;
    }
    catch (...)
    {
        std::cout<<">> 未知的异常! "<<std::endl;
    }

    std::cout<<std::endl;
    std::cout<<"...end..."<<std::endl;
    char c; std::cin>>c;
    
    return 0; 

}


//计算圆周率PI
TLargeFloat GetBorweinPI(unsigned long uiDigitsLength)
{
    //Borwein四次迭代式：
    //    初值： a0=6-4*2^0.5; y0=2^0.5-1;    
    //    重复计算： y(n+1)=[1-(1-y^4)^0.25]/[1+(1-y^4)^0.25];
    //               y=y(n+1);
    //               a(n+1)=a*(1+y)^4-2^(2*n+3)*y*(1+y+y*y);
    //    最后计算： PI=1/a;

    //迭代次数和10进制有效精度  每迭代一次有效精度增长4倍
    //0 1  2  3   4    5    6   ...
    //0 8 40 170 693 2789 11171 ...
    long loop_count; 
    if (uiDigitsLength<=8)
        loop_count=1;
    else if (uiDigitsLength<=40)
        loop_count=2;
    else if (uiDigitsLength<=170)
        loop_count=3;
    else if (uiDigitsLength<=693)
        loop_count=4;
    else if (uiDigitsLength<=2789)
        loop_count=5;
    else 
    {
        loop_count=6;
        long DigitsLength=11171;
        while (DigitsLength<(long)uiDigitsLength)
        {
            DigitsLength*=4;//四阶收敛
            ++loop_count;
        }
    }
    std::cout<<"        ... 开始计算PI,Borwein四次迭代式, 总循环次数："<<loop_count<<"  ...   "<<std::endl;

    TLargeFloat::TDigits sCount(uiDigitsLength);//计算采用的十进制精度

    TLargeFloat a(0.0,sCount),y(0.0,sCount);
    TLargeFloat temp(0.0,sCount);
    TLargeFloat pow(2*2*2,sCount); //2^(2*n+3); (n=0);

    //a0=6-4*2^0.5; 
    temp=2;
    temp.Sqrt();
    a=6-4*temp;
    //y0=2^0.5-1;    
    y=temp-1;


    TLargeFloat tmpa(0.0,sCount);
    for (long i=0;i<loop_count;++i)
    {
        std::cout<<"        ... 计算PI "<<uiDigitsLength<<" 位的当前循环次数："<<i+1<<"/"<<loop_count<<"  ...   "<<std::endl;
       
        temp=1-sqr(sqr(y));
        temp=revsqrt(revsqrt(temp));//等价于 temp=sqrt(sqrt(temp));
        y=(1-temp)/(1+temp);
        temp=sqr(1+y);
        a=a*sqr(temp)-pow*y*(temp-y);

        pow*=4;
    }
    std::cout<<"        ... 计算PI的循环结束，准备输出  ...   "<<std::endl; 
    return 1/a;//PI=1/a;
}


TLargeFloat GetAGMPI(unsigned long uiDigitsLength)
{
    //AGM二次迭代式:
    //初值:a=x=1 b=1/sqrt(2) c=1/4
    //重复计算:y=a a=(a+b)/2 b=sqrt(by) c=c-x(a-y)^2 x=2x
    //最后:pi=(a+b)^2/(4c)

    //迭代次数和10进制有效精度  每迭代一次有效精度增长2倍
    //0  1  2  3   4    5    6    7     8  ...
    //1 39 82 169 344  693 1391 2788  5582 ...
    long loop_count; 
    if (uiDigitsLength<=39)
        loop_count=1;
    else if (uiDigitsLength<=82)
        loop_count=2;
    else if (uiDigitsLength<=169)
        loop_count=3;
    else if (uiDigitsLength<=344)
        loop_count=4;
    else if (uiDigitsLength<=693)
        loop_count=5;
    else if (uiDigitsLength<=1391)
        loop_count=6;
    else if (uiDigitsLength<=2788)
        loop_count=7;
    else 
    {
        loop_count=8;
        long DigitsLength=5582;
        while (DigitsLength<(long)uiDigitsLength)
        {
            DigitsLength*=2;//2阶收敛
            ++loop_count;
        }
    }
    std::cout<<"        ... 开始计算PI,AGM二次迭代式, 总循环次数："<<loop_count<<"  ...   "<<std::endl;

    TLargeFloat::TDigits sCount(uiDigitsLength);//计算采用的十进制精度

    TLargeFloat c(8,sCount);  c.RevSqrt(); c.Chs();  //c=-sqrt(1/8);
    TLargeFloat a(c); a*=(-3); a+=1;                 //a=1-3*c;
    TLargeFloat b(a); b.Sqrt();                      //b=sqrt(a);
    TLargeFloat e("-0.625"); e += b;                 //e=b-0.625
    b *= 2;
    c += e;
    a += e;
    long npow=4;
    for (long i=0;i<loop_count;++i)
    {
        std::cout<<"        ... 计算PI "<<uiDigitsLength<<" 位的当前循环次数："<<i+1<<"/"<<loop_count<<"  ...   "<<std::endl;
        e  = a;
        e += b;
        e /= 2;
        b *=a;
        b.Sqrt();
        e -= b;
        b *= 2;
        c -= e;
        a  = e; 
        a += b;
        npow *= 2;
    }
    std::cout<<"        ... 计算PI的循环结束，准备输出  ...   "<<std::endl; 
    e *= e;
    e /= 4;
    a += b;

    c *=a;
    c -=e;
    c *=npow;
    
    a *=a;
    a-=e;
    e/=2;
    a-=e;

    a/=c;

    return  a;
}
//
///////////////////////////////////////////////////////////////////////////////////