tlargefloat.cpp

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

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// TLargeFloat.cpp: implementation of the TLargeFloat class.
//   超高精度浮点数类TLargeFloat
//////////////////////////////////////////////////////////////////////

#include "TLargeFloat.h"
#include <algorithm>
//#include "assert.h"
#include <math.h>

//#define MyDebugAssert(b) { if (!(b)) __asm{ int 3 } }

//用以控制运算精度的工具
class TDigitsLengthCtrl
{
private:
    std::vector<TLargeFloat*> m_list;
public:
    inline void add_to_ctrl(TLargeFloat* var) { m_list.push_back(var); }
    void SetDigitsLength(const long uiDigitsLength) 
    { 
        long size=m_list.size();
        for (long i=0;i<size;++i)
            m_list[i]->SetDigitsLength(uiDigitsLength);
    }
};

//自动精度控制
class TDigitsLengthAutoCtrl:public TDigitsLengthCtrl
{
private:
    std::vector<long>    m_DigitsLengthList;
    long   step;
public:
    TDigitsLengthAutoCtrl(long DigitsMul,long DigitsLength0,long MaxDigitsLength)//(收敛系数,初始精度,最终需要精度)
    {
        //MyDebugAssert(DigitsLength0>0);
        //得到最佳的精度梯度
        long curDigitsLength=MaxDigitsLength;
        while (curDigitsLength>DigitsLength0)
        {
            m_DigitsLengthList.push_back(curDigitsLength);

            if ((DigitsLength0<=1)&&(curDigitsLength==1)) break;
            curDigitsLength=(curDigitsLength+DigitsMul-1)/DigitsMul;
        }
        step=get_runcount()-1;
    }
    inline long get_runcount() { return m_DigitsLengthList.size(); }
    inline void SetNextDigits()
    {
        //MyDebugAssert(step>=0);
        long DigitsLength=m_DigitsLengthList[step];
        if (step>0)
            DigitsLength+=4;
        SetDigitsLength(DigitsLength);
        --step;
        if (step<0) step=0;
    }
};

/////////////////////////////////////////////////////////////////////////////////////////////

void TLargeFloat_CtrlExponent(TLargeFloat::TArray& v,TMaxInt iMoveRightCount)
{
    //iMoveCountBig 有可能存不下 但现在的体系下因右移超过值域比较难
    long iMoveCountBig=(long)(iMoveRightCount/TLargeFloat::em10Power);
    long iMoveCountSmall=(long)(iMoveRightCount%TLargeFloat::em10Power);

    //大的右移
    if (iMoveCountBig>0)
    {
        for (long i=v.size()-1;i>=iMoveCountBig;--i)
            v[i]=v[i-iMoveCountBig];
        for (long j=0;j<iMoveCountBig;++j)
            v[j]=0;
    }
    //小的右移
    if (iMoveCountSmall>0)
    {
        TInt32bit iDiv=1;
        for (long j=0;j<iMoveCountSmall;++j) iDiv*=10;

        long v_size=v.size();
        for (long i=iMoveCountBig;i<v_size-1;++i)
        {
            v[i+1]+=(v[i]%iDiv)*TLargeFloat::emBase;
            v[i]/=iDiv;
        }
        if (v_size>=1)
            v[v_size-1]/=iDiv;
    }
    
}

    void TLargeFloat_ArrayAddTestLast(TInt32bit* x,long xsize)
    {
        for (long i=xsize-1;i>0;--i)
        {
            if (x[i]>=TLargeFloat::emBase)//进位
            {
                ++x[i-1];
                x[i]-=TLargeFloat::emBase;
            }
            else
                break;
        }   
    }

void TLargeFloat_ArrayAdd(TInt32bit* result,long rsize,const TInt32bit* y,long ysize)
{
    //MyDebugAssert(rsize>=ysize);
    //MyDebugAssert(ysize>0);
    TInt32bit* ri=&result[rsize-ysize];
    long i;
    for (i=ysize-1;i>0;--i)
    {
        TInt32bit r=ri[i]+y[i];
        if (r>=TLargeFloat::emBase)//进位
        {
            ++ri[i-1];
            r-=TLargeFloat::emBase;
        }
        ri[i]=r;
    }    
    ri[0]+=y[0];
  
    TLargeFloat_ArrayAddTestLast(result,rsize-ysize+1);
}

void TLargeFloat_ArraySub(TInt32bit* result,long rsize,const TInt32bit* y,long ysize)
{
    //MyDebugAssert(rsize>=ysize);
    //MyDebugAssert(ysize>0);
    TInt32bit* ri=&result[rsize-ysize];
    long i;
    for (i=ysize-1;i>0;--i)
    {
        TInt32bit r=ri[i]-y[i];
        if (r<0)//借位
        {
            --ri[i-1];
            r+=TLargeFloat::emBase;
        }
        ri[i]=r;
    }    
    ri[0]-=y[0];
  
    for (i=rsize-ysize;i>0;--i)
    {
        if (result[i]<0)//借位
        {
            --result[i-1];
            result[i]+=TLargeFloat::emBase;
        }
        else
            break;
    }   
}


        TInt32bit _inti_csMaxMuliValue()
        {
            double rb=1; double mb=1;
            for (long i=0;i<100/TLargeFloat::em10Power;++i)
            {
                mb*=(1.0/TLargeFloat::emBase);
                rb+=mb;
            }
            return (TInt32bit)( TInt32bit_MAX_VALUE*(1.0/(TLargeFloat::emBase-1))/rb ) -1;
        }
    //不超界所允许的乘数值M满足: (emBase-1)*M + (emBase-1)*M/B + (emBase-1)*M/B^2 + (emBase-1)*M/B^3 ... <=TInt32bit_MAX_VALUE
    static const TInt32bit csMaxMuliValue=_inti_csMaxMuliValue();
    static const TInt32bit csMaxMulAddValue = csMaxMuliValue - (TLargeFloat::emBase-1);

void TLargeFloat_ArrayMUL_ExE(TInt32bit* result,long rminsize,const TInt32bit* x,long xsize,const TInt32bit* y,long ysize)
{
    //N*N复杂度的简单乘法实现
    for (long k=0;k<rminsize;++k) result[k]=0;
    while ((ysize>0)&&(y[0]==0)) { --ysize; ++y; ++result; --rminsize;}
    //while ((xsize>0)&&(x[0]==0)) { --xsize; ++x; ++result; --rminsize;}
    while ((ysize>0)&&(y[ysize-1]==0)) { --ysize; }
    //while ((xsize>0)&&(x[xsize-1]==0)) { --xsize; }

    long add_sum=0;
    for (long i=0;i<xsize;++i)
    {
        long xi=x[i];
        if (xi>0)
        {
            TInt32bit* resulti=&result[i+1];
            long y_limit=std::min(rminsize-1-i,ysize);
            for (long j=0;j<y_limit;++j)
            {
                resulti[j]+=(xi*y[j]);
            }
            add_sum+=xi;
            if (add_sum>csMaxMulAddValue)
            {
                for (long k=i+y_limit;k>0;--k)
                {
                    TInt32bit v=result[k];
                    if (v>=TLargeFloat::emBase)
                    {
                        result[k-1]+=v/TLargeFloat::emBase;
                        result[k]= v%TLargeFloat::emBase;
                    }
                }
                add_sum=0;
            }
        }
    }
    if (add_sum>0)
    {
        for (long k=rminsize-1;k>0;--k)
        {
            TInt32bit v=result[k];
            if (v>=TLargeFloat::emBase)
            {
                result[k-1]+=v/TLargeFloat::emBase;
                result[k]= v%TLargeFloat::emBase;
            }
        }
    }
}

//数组乘法  核心
void ArrayMUL(TInt32bit* result,long rminsize,const TInt32bit* x,long xsize,const TInt32bit* y,long ysize);


void ArrayMUL_Dichotomy(TInt32bit* result,long rminsize,const TInt32bit* x,const TInt32bit* y,long MulSize)
{

    //二分法的乘法实现
    //   x*y=(a*E+b)*(c*E+d)=(E*E-E)*ac + E*(a+b)*(c+d) -(E-1)*bd

    //   temp_a_add_b=a+b;
    //   temp_c_add_d=c+d;
    //   result=E*temp_a_add_b*temp_c_add_d;

    //   temp_ac=a*c;
    //   result+=E*E*temp_ac;
    //   result-=E*temp_ac;

    //   temp_bd=b*d;
    //   result+=temp_bd;
    //   result-=E*temp_bd;
    // 
    long i;

    const long E=((MulSize+1)>>1);
    //if (tmpData+2*(E+1)>end_tmpData)  __asm int 3
    const long acSize=MulSize-E;
    const long acErrorSize=E-acSize; 

    const TInt32bit* a=x;
    const TInt32bit* b=&x[acSize];
    const TInt32bit* c=y;
    const TInt32bit* d=&y[acSize];

    TLargeFloat::TArray _tmpData((E+1)*2);
    TInt32bit* tmpData=&_tmpData[0];

    TInt32bit* temp_a_add_b=&tmpData[0];
    TInt32bit* temp_c_add_d=&tmpData[(E+1)];

    temp_a_add_b[0]=0;     
    for (i=0;i<E;++i) temp_a_add_b[i+1]=b[i];
    TLargeFloat_ArrayAdd(temp_a_add_b,1+E,a,E-acErrorSize);
    temp_c_add_d[0]=0;
    for (i=0;i<E;++i) temp_c_add_d[i+1]=d[i];
    TLargeFloat_ArrayAdd(temp_c_add_d,1+E,c,E-acErrorSize);

    const long abcdPos=(2*MulSize-E-2*(E+1));
    for (i=0;i<abcdPos;++i) result[i]=0;
    ArrayMUL(&result[abcdPos],std::min(rminsize-abcdPos,2*(E+1)),temp_a_add_b,E+1,temp_c_add_d,E+1);
    for (i=abcdPos+2*(E+1);i<rminsize;++i) result[i]=0;

    TInt32bit* temp_ac=&tmpData[0];
    const long acminsize=std::min(rminsize,acSize*2);
    ArrayMUL(temp_ac,acminsize,a,acSize,c,acSize);
    TLargeFloat_ArrayAdd(result,acminsize,temp_ac,acminsize);
    long sub_ac_size=acSize*2;
    if (E+sub_ac_size>rminsize) sub_ac_size=rminsize-E;
    TLargeFloat_ArraySub(result,E+sub_ac_size,temp_ac,sub_ac_size);

    const long add_bdPos=2*MulSize-2*E;
    long sub_bd_size=2*E;
    const long sub_bdPos=add_bdPos-E;
    if (sub_bdPos+sub_bd_size>rminsize) sub_bd_size=rminsize-sub_bdPos;
    if (sub_bd_size>0)
    {
        TInt32bit* temp_bd=&tmpData[0];
        ArrayMUL(temp_bd,sub_bd_size,b,E,d,E);
        long add_bd_size=2*E;
        if (add_bdPos+add_bd_size>rminsize) add_bd_size=rminsize-add_bdPos;
        if (add_bd_size>0)
          TLargeFloat_ArrayAdd(result,add_bdPos+add_bd_size,temp_bd,add_bd_size);
        TLargeFloat_ArraySub(result,sub_bdPos+sub_bd_size,temp_bd,sub_bd_size); 
    }
}


void ArrayMUL_DichotomyPart(TInt32bit* result,long rminsize,const TInt32bit* x,long xsize,const TInt32bit* y,long ysize)
{
    //二分法的辅助函数
    ///////////////////////////////////////////
    //x*y=(a*E+b)*y=a*y*E+b*y;

    if (xsize==ysize)
    {
        ArrayMUL_Dichotomy(result,rminsize,x,y,ysize);
        return ;
    }

    if (xsize<ysize)
    {
        std::swap(xsize,ysize);
        std::swap(x,y);
    }

    const long E=ysize;
    const long asize=xsize-E;
    const TInt32bit* a=&x[0];
    const TInt32bit* b=&x[asize];

 
    ArrayMUL(result,rminsize,a,asize,y,ysize);
    for (long i=asize+ysize;i<rminsize;++i) result[i]=0;
    const long byPos=asize;
    if (byPos<rminsize)
    {
        const long byminsize=rminsize-byPos;
        TLargeFloat::TArray tmpData(byminsize);
        TInt32bit* tmp_by=&tmpData[0];
        ArrayMUL_Dichotomy(tmp_by,byminsize,b,y,E);
        TLargeFloat_ArrayAdd(result,rminsize,tmp_by,byminsize);
    }
}


/////////////////////////////////////////////////////////////////////////////////////////////

const long csTMaxInt_Digits=(long)(log10((long double)TMaxInt_MAX_VALUE)+1);

TLargeFloat::TLargeFloat(const SelfType& Value)
:m_Digits(Value.m_Digits)
{
    m_Exponent=Value.m_Exponent;
    m_Sign=Value.m_Sign;
}
    
TLargeFloat::TLargeFloat(const long double DefultValue)
:m_Digits(TDigits(emLongDoubleDigits+1).GetDigitsArraySize(),0)
{
    fToLargeFloat(DefultValue);
}

TLargeFloat::TLargeFloat(const long double DefultValue,const TDigits& DigitsLength)
:m_Digits(DigitsLength.GetDigitsArraySize(),0)
{
    fToLargeFloat(DefultValue);
}

TLargeFloat::TLargeFloat(const char* strValue)
:m_Digits(TDigits(csTMaxInt_Digits).GetDigitsArraySize(),0)
{
    sToLargeFloat(std::string(strValue));
}
TLargeFloat::TLargeFloat(const char* strValue,const TDigits& DigitsLength)
:m_Digits(DigitsLength.GetDigitsArraySize(),0)
{
    sToLargeFloat(std::string(strValue));
}
TLargeFloat::TLargeFloat(const std::string& strValue)
:m_Digits(TDigits(csTMaxInt_Digits).GetDigitsArraySize(),0)
{
    sToLargeFloat(strValue);
}
TLargeFloat::TLargeFloat(const std::string& strValue,const TDigits& DigitsLength)
:m_Digits(DigitsLength.GetDigitsArraySize(),0)
{
    sToLargeFloat(strValue);
}


void TLargeFloat::Zero()
{
    if (m_Sign==0) return;
    m_Sign=0;
    m_Exponent=0;
    long old_size=m_Digits.size();
    for (long i=0;i<old_size;++i)
        m_Digits[i]=0;
}

void TLargeFloat::Canonicity()
{
    //规格化 小数部分
    //规格化前允许:m_Digits[0]>=emBase 也就是小数部分的值大于等于1.0;
    //规格化前允许:m_Digits[]前面有多个0值 也就是小数部分的值小于0.1;
    //规格化后的数满足的条件：小数部分的值 小于1.0，大于等于0.1;

    if (m_Sign==0)
    {
        return;
    }
    long old_size=m_Digits.size();
    if (old_size<=0) return; //不可能发生:)

    if ( (m_Digits[0]<emBase)&&(m_Digits[0]>=(emBase/10)) ) return;

    //处理右移
    while (m_Digits[0]>=emBase)//是否需要右移1
    {
        m_Exponent+=em10Power;
        for (long i=old_size-1;i>=2;--i)
            m_Digits[i]=m_Digits[i-1];
        if (old_size>=2) m_Digits[1]=m_Digits[0] % emBase;
        m_Digits[0]/=emBase;
    }

    ////////////////////////////////////
    //考虑左移

    //大的左移
    long iMoveCountBig=old_size;
    for (long i=0;i<old_size;++i)
    {
        if (m_Digits[i]!=0)
        {
            iMoveCountBig=i;
            break;
        }
    }
    if (iMoveCountBig==old_size)
    {
        //as Zero();
        m_Sign=0;
        m_Exponent=0;
        return;
    }

    if (iMoveCountBig>0)
    {
        m_Exponent-=(iMoveCountBig*em10Power);
        for (long j=0;j<old_size-iMoveCountBig;++j)
            m_Digits[j]=m_Digits[j+iMoveCountBig];
        for (long k=old_size-iMoveCountBig;k<old_size;++k)
            m_Digits[k]=0;
    }

    //小的左移
    TInt32bit iMoveMul=1;
    long iMoveCountSmall=0;
    while (iMoveMul*m_Digits[0]<(emBase/10))
    {
        iMoveMul*=10;
        ++iMoveCountSmall;
    }
    if (iMoveCountSmall>0) 
    {
        m_Exponent-=iMoveCountSmall;
        for (long i=old_size-1;i>=0;--i)
            m_Digits[i]*=iMoveMul;
        for (long j=old_size-1;j>0;--j)
        {
            TInt32bit value=m_Digits[j];
            m_Digits[j]=value%emBase;
            m_Digits[j-1]+=value/emBase;
        }
    }
}

void TLargeFloat::iToLargeFloat(TMaxInt iValue)
{
    Zero();
    long new_size=TDigits(csTMaxInt_Digits).GetDigitsArraySize();
    if ((long)m_Digits.size()<new_size)
        m_Digits.resize(new_size,0);
    if (iValue==0)
        return;
    else if (iValue>0)
        m_Sign=1;
    else
    {
        m_Sign=-1;
        iValue=-iValue;
    }

    //无误差转换
    TMaxInt tmp_iValue=iValue;
    long n=0;
    for (;tmp_iValue!=0;++n)
        tmp_iValue/=emBase;
    m_Exponent=n*em10Power;
    for (long i=0;i<n;++i)
    {
        m_Digits[n-1-i]=(TInt32bit)(iValue%emBase);
        iValue/=emBase;
    }
    Canonicity();
}

void TLargeFloat::fToLargeFloat(long double fValue)
{
    TMaxInt iValue=(TMaxInt)fValue;
    if (iValue==fValue)
    {
        iToLargeFloat(iValue);
        return;
    }