tlargefloat.cpp

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

        TMaxInt iDiv=iValue;
        long size=m_Digits.size();
        TMaxInt cur_value=m_Digits[0];
        for (long i=0;i<size-1;++i)
        {
            m_Digits[i]=(TInt32bit)(cur_value/iDiv);
            cur_value=m_Digits[i+1]+(cur_value%iDiv)*emBase;
        }
        if (size>=1)
            m_Digits[size-1]=(TInt32bit)(cur_value/iDiv);
    }

    Canonicity();
}

TLargeFloat::SelfType& TLargeFloat::operator /= (long double  fValue)
{
    TMaxInt iValue=(TMaxInt)fValue;
    if (iValue==fValue)
    {
        DivInt(iValue);
        return *this;
    }
    else
        return (*this)/=TLargeFloat(fValue);
}

TLargeFloat::SelfType& TLargeFloat::operator /= (const SelfType& Value)
{
    SelfType x(Value);
    x.Rev();//x=1/x;
    return (*this)*=x;
}

void TLargeFloat::Rev()//求倒数
{
    //求1/a,
    // 则有a=1/x;  y=a-1/x;
    // 求导得y'=1/x^2;
    // 代入牛顿方程 x(n+1)=x(n)-y(x(n))/y'(x(n));
    // 有迭代式 x_next=(2-a*x)*x; //可证明：该公式为2阶收敛公式
    // 证明:令x=1/a+dx;
    // 则有x_next-1/a=(2-a*(1/a+dx))*(1/a+dx)-1/a
    //               =(-a)*dx^2
    // 证毕.

    if (this->m_Sign==0) throw TException("ERROR:TLargeFloat::Rev()");
    TExpInt oldExponent=m_Exponent;

    SelfType a(*this);
    a.m_Exponent=0;
    SelfType x(0.0,TLargeFloat::TDigits(a.GetDigitsLength()));
    x=1.0/a.AsFloat();//初始迭代值

    TDigitsLengthAutoCtrl dlCtrl(2,emLongDoubleDigits-1,(a.m_Digits.size()*TLargeFloat::em10Power));
    SelfType ta(0.0,TLargeFloat::TDigits(a.GetDigitsLength()));
    SelfType temp(0.0,TLargeFloat::TDigits(a.GetDigitsLength()));
    dlCtrl.add_to_ctrl(&ta);
    dlCtrl.add_to_ctrl(&x);
    dlCtrl.add_to_ctrl(&temp);
    for (long i=0;i<dlCtrl.get_runcount();++i)
    {
        ta=a;
        dlCtrl.SetNextDigits();
        //x=x*2-x*x*ta;
        temp=x;
        temp*=temp;
        temp*=ta;
        x*=2;
        x-=temp;
    }

    x.m_Exponent-=oldExponent;
    this->Swap(x);
}


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

void TLargeFloat::RevSqrt()//
{
    //求1/a^0.5,
    // 则有a=1/x^2;  y=a-1/x^2;
    // 求导得y'=2/x^3;
    // 代入牛顿方程 x(n+1)=x(n)-y(x(n))/y'(x(n));
    // 有迭代式 x_next=(3-a*x*x)*x/2; //可证明：该公式为2阶收敛公式
    // 证明:令x=1/a^0.5+dx;
    // 则有x_next-1/a^0.5=(3-a*(1/a^0.5+dx)*(1/a^0.5+dx))*(1/a^0.5+dx)/2 - 1/a^0.5
    //                   =-(1.5/a^0.5)*dx^2-0.5*a*dx^3;
    // 证毕.
    if (this->m_Sign<0) throw TException("ERROR:TLargeFloat::RevSqrt()");
    if (this->m_Sign==0) return;

    SelfType a(*this);
    TMaxInt sqrExponent=a.m_Exponent.AsInt()/2;
    a.m_Exponent-=(sqrExponent*2);


    SelfType x(0.0,TDigits(a.GetDigitsLength()));
    x=1.0/sqrt(a.AsFloat());//初始迭代值

    TDigitsLengthAutoCtrl dlCtrl(2,emLongDoubleDigits-1,(a.m_Digits.size()*TLargeFloat::em10Power));
    SelfType ta(0.0,TLargeFloat::TDigits(a.GetDigitsLength()));
    SelfType temp(0.0,TLargeFloat::TDigits(a.GetDigitsLength()));
    dlCtrl.add_to_ctrl(&ta);
    dlCtrl.add_to_ctrl(&x);
    dlCtrl.add_to_ctrl(&temp);
    for (long i=0;i<dlCtrl.get_runcount();++i)
    {
        ta=a;
        dlCtrl.SetNextDigits();
        //x=(3-x*x*ta)*x/2;
        temp=x;
        temp*=temp;
        temp*=ta;
        temp.Chs();
        temp+=3;
        x*=temp;
        x/=2;
    }

    x.m_Exponent-=sqrExponent;
    this->Swap(x);
}

////////////////////////////
//*

void TLargeFloat::MulInt(TMaxInt iValue)
{
    if (this->m_Sign==0) return;
    if (iValue==0) 
    {
        this->Zero();
        return;
    }
    else if (iValue<0)
    {
        iValue=-iValue;
        this->Chs();
        //continue
    }

    if (iValue>csMaxMuliValue)
    {
        TLargeFloat y(0.0);
        y.iToLargeFloat(iValue);
        (*this)*=y;
        return;
    }

    //乘以一个整数
    TInt32bit Value=(TInt32bit)iValue;
    long size=m_Digits.size();
    for (long i=size-1;i>=0;--i)
    {
        m_Digits[i]*=Value;
    }
    for (long j=size-1;j>=1;--j)
    {
        if (m_Digits[j]>=emBase)
        {
            m_Digits[j-1]+=m_Digits[j]/emBase;
            m_Digits[j]%=emBase;
        }
    }
    Canonicity();
}
TLargeFloat::SelfType& TLargeFloat::operator *= (long double  fValue)
{
    TMaxInt iValue=(TMaxInt)fValue;
    if (iValue==fValue)
    {
        MulInt(iValue);
        return *this;
    }
    else
    {
        TLargeFloat temp(fValue);
        return (*this)*=temp;
    }
}

inline void getBestMulSize(const TInt32bit*& x,long& xsize,const TInt32bit*& y,long& ysize)
{
    long bxsize=xsize;
    long bysize=ysize;

    long i;
    for (i=xsize-1;i>=0;--i)
    {
        if (x[i]==0) 
            --bxsize;
        else
            break;
    }
    for (i=ysize-1;i>=0;--i)
    {
        if (y[i]==0) 
            --bysize;
        else
            break;
    }

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

    const long csABigMulSize=1000;
    if ( (bysize<csABigMulSize)||(bxsize*2>=bysize*3)||(bxsize>ysize) )
        ysize=bysize;
    else 
        ysize=bxsize;
}

TLargeFloat::SelfType& TLargeFloat::operator *= (const SelfType& Value)
{
    long xsize=m_Digits.size();
    long ysize=Value.m_Digits.size();
    if ((m_Sign==0)||(Value.m_Sign==0)) 
    {
        if (m_Sign!=0)
            this->Zero();
        if (xsize<ysize)
            m_Digits.resize(ysize,0);
        return *this;
    }
    if (xsize<ysize)
        m_Digits.resize(ysize,0);
    
    long rminsize=std::max(xsize,ysize)+2;//乘法结果需要的位数

    const TInt32bit* x=&m_Digits[0];
    const TInt32bit* y=&Value.m_Digits[0];
    getBestMulSize(x,xsize,y,ysize);
    if (rminsize>xsize+ysize)
        rminsize=xsize+ysize;

    TArray tempArray(rminsize);
    ArrayMUL(&tempArray[0],rminsize,x,xsize,y,ysize);

    xsize=m_Digits.size();
    long resultsize=std::min(xsize,rminsize);
    long i;
    for (i=0;i<resultsize;++i)  this->m_Digits[i]=tempArray[i];
    for (i=resultsize;i<xsize;++i) this->m_Digits[i]=0;
    m_Exponent+=Value.m_Exponent;
    m_Sign*=Value.m_Sign;
    Canonicity();
    return *this;
}


////////////////////////////////////////////////////////////
//快速傅立叶变换实现大整数乘法 
//  FFT Mul

struct TComplex //复数类型    
{
    double  i;
    double  r;
};

//初始化

void FFTInitial(TComplex* w,long* lst,const long n)
{
    //计算W^i , i in [0--n)
    const double PI=3.1415926535897932385;
    const double seta=2.0*PI/n;
    long i,j;

    long n2=n>>2;
    for (i=0;i<=n2;++i)
    {
        w[i].i=sin(seta*i);
        w[i].r=cos(seta*i);
    }
    for (j=(n2+1),i=(n2-1); i>=0 ;--i,++j)
    {
        w[j].i= w[i].i;
        w[j].r=-w[i].r;
    }

    //计算排序用数组
    lst[0]=0; 
    i=1;
    long m;
    for (n2=n>>1,m=2; n2>0 ;n2=n2>>1,m =m<<1) 
    {
        for (long j=0;i<m;++i,++j)
        {
            lst[i]=lst[j]+n2;
        }
    }
}

//快速傅立叶变换
void FFT(TComplex* a,const TComplex* w,const long n)
{
    long n2=1;  
    long n1=n >> 1;
    for (long per=1;per<n;)
    {
        long m=n2;  long j=0;  long k=n2;
        per<<=1;
        while (k<n) 
        {
            for (long i2=0; j<m ;++j,++k,i2+=n1)
            {
                TComplex tmpb=a[j];
                double   tmp1=a[k].r*w[i2].r-a[k].i*w[i2].i;
                double   tmp2=a[k].r*w[i2].i+a[k].i*w[i2].r;
                a[j].r=tmpb.r+tmp1;
                a[j].i=tmpb.i+tmp2;
                a[k].r=tmpb.r-tmp1;
                a[k].i=tmpb.i-tmp2;
            }
            m+=per;  
            j+=n2;  
            k+=n2;
        }
        n2=per;
        n1>>=1;
    }
}


//反向快速傅立叶变换
void InvFFT(TComplex* a,const TComplex* w,const long n)
{
    for (long i=0;i<n;++i)
        a[i].i=-a[i].i;
    FFT(a,w,n);
}

void FFT_Convolution(const TComplex* a,TComplex* b,const long n)
{
    for (long i=0;i<n;++i)
    {
        double br=a[i].r*b[i].r - a[i].i*b[i].i;
        b[i].i=a[i].r*b[i].i + a[i].i*b[i].r;
        b[i].r=br;
    }
}

    void FFT_IntArrayCopyToComplexArray(const TInt32bit* ACoef,const long ASize,TComplex* result,const long n,const long* lst)
    {
        for (long i=0;i<n;++i)
        {
            long ai=lst[i];
            if (ai<ASize)
                result[i].r=ACoef[ai];  
            else
                result[i].r=0.0;  
            result[i].i=0.0;
        }
    }

    void FFT_ComplexArrayToIntArray(const TComplex* cmx,const long n,TInt32bit* result,const long rMinSize)
    {
        const double tmpRevN=1.0/n;
        double FFT_int_base;
        TLargeFloat::TArray tempResult;
        TInt32bit* dst=0;
        long dstsize=0;
        {   
            FFT_int_base=(long)TLargeFloat::emBase;
            dstsize=rMinSize;
            dst=result;
        }
        const double tmpRevFFT_int_base=1.0/FFT_int_base;

        double MaxWorstEr=0;
        double tmp=0.0;
        for (long i=std::min(n,dstsize)-1;i>=1;--i)
        {
            double num =cmx[i-1].r*tmpRevN+tmp;
            tmp=floor((num+0.499)*tmpRevFFT_int_base);
            dst[i]=(TInt32bit)(num+0.499-tmp*FFT_int_base);

            double newWorstEr=abs(num-tmp*FFT_int_base-dst[i]);//记录最大误差
           if (newWorstEr>MaxWorstEr)
                MaxWorstEr=newWorstEr;  
        }
        if (dstsize>0)
            dst[0]=(TInt32bit)(tmp+0.1);
        if (MaxWorstEr>=0.499) 
            throw TLargeFloat::TException("ERROR:FFT's FFT_ComplexArrayToIntArray(); "); //误差过大
    }


//FFT平方
void SquareFFT(const TInt32bit* ACoef,const long ASize,TInt32bit* CCoef,const long CMinSize)
{
    long ansize=ASize+ASize;
    long n=1;
    while (n<ansize) n<<=1;

    std::vector<TComplex>   _w((n>>1)+1);   
    std::vector<long>       _lst(n);        
    TComplex* w=&_w[0];
    long*     lst=&_lst[0];
    FFTInitial(w,lst,n);

    std::vector<TComplex>   _a(n);  
    TComplex* a=&_a[0];
    FFT_IntArrayCopyToComplexArray(ACoef,ASize,a,n,lst);

    FFT(a,w, n);
    FFT_Convolution(a,a,n);

    std::vector<TComplex>   _tmpa(n);
    _tmpa.swap(_a);
    TComplex* tmpa=&_tmpa[0]; a=&_a[0];
    for (long i=0;i<n;++i)  a[i]=tmpa[lst[i]];  
    _tmpa.clear();
    _lst.clear();
    InvFFT(a,w,n);
    _w.clear();
    FFT_ComplexArrayToIntArray(a,n,CCoef,CMinSize);
 }


//FFT乘法
void MulFFT(const TInt32bit* ACoef,const long ASize,const TInt32bit* BCoef,const long BSize,TInt32bit* CCoef,const long CMinSize)
{
    long absize=(ASize+BSize);
    long n=1;
    while (n<absize) n<<=1;

    std::vector<TComplex>   _w((n>>1)+1);
    std::vector<long>       _lst(n);
    TComplex* w=&_w[0];
    long*     lst=&_lst[0];
    FFTInitial(w,lst,n);

    std::vector<TComplex>   _a(n);
    TComplex* a=&_a[0];
    FFT_IntArrayCopyToComplexArray(ACoef,ASize,a,n,lst);
    std::vector<TComplex>   _b(n);
    TComplex* b=&_b[0];
    FFT_IntArrayCopyToComplexArray(BCoef,BSize,b,n,lst);

    FFT(b,w, n);
    FFT(a,w, n);
    FFT_Convolution(a,b,n);
    for (long i=0;i<n;++i)  a[i]=b[lst[i]];  
    _b.clear();
    _lst.clear();
    InvFFT(a,w,n);
    _w.clear();
    FFT_ComplexArrayToIntArray(a,n,CCoef,CMinSize);
}


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

const long csBestMUL_ExE_size           =1024       /TLargeFloat::em10Power;
const long csBestMUL_Dichotomy_size     =4096       /TLargeFloat::em10Power;
const long csBestMul_FFT4Max_size       =16777216   /TLargeFloat::em10Power;
void ArrayMUL(TInt32bit* result,long rminsize,const TInt32bit* x,long xsize,const TInt32bit* y,long ysize)
{
    if (xsize<ysize)
    {
        std::swap(xsize,ysize);
        std::swap(x,y);
    }
    if (rminsize<ysize)
    {
        xsize=rminsize;
        ysize=rminsize;
    }
    if (ysize<=csBestMUL_ExE_size)  //
        TLargeFloat_ArrayMUL_ExE(result,rminsize,x,xsize,y,ysize);
    else  if (ysize<=csBestMUL_Dichotomy_size)
        ArrayMUL_DichotomyPart(result,rminsize,x,xsize,y,ysize);
    else if ((xsize+ysize)<=(csBestMul_FFT4Max_size*2))
    {
        if (x==y)
            SquareFFT(x,xsize,result,rminsize);
        else
            MulFFT(x,xsize,y,ysize,result,rminsize);
    }
    else
        ArrayMUL_DichotomyPart(result,rminsize,x,xsize,y,ysize);
}



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

void LargeFloat_UnitTest()
{
    //单元测试

    //测试先行... 但作者太懒...
    //todo:

    TLargeFloat x(200,TLargeFloat::TDigits(100));
    TLargeFloat y(12345.678987654321,TLargeFloat::TDigits(100));
    //TLargeFloat z("-12345.678987654321e-100",TLargeFloat::TDigits(100));

    TLargeFloat z("99999999999999999999999999999999999999999999",TLargeFloat::TDigits(200));
//
}