高精度浮点数.c

一个通用的高精度浮点数算法的实现,可以学习一下.

    if (m_Sign==0) return 0;
    long double result=m_Sign;
    result*=pow(10,double(m_Exponent.AsInt()));
    
    long double r=1;
    long double Sum=0;
 int m_CalcSize=m_Digits.size();
    for (int i=0;i<m_CalcSize;++i)//得到小数位
    {
        r/=emBase;
        if (r==0) break; 
        Sum+=(m_Digits[i]*r);
    }       
    return result*Sum;
}


void TLargeFloat::MoveLeft10Power(TMaxInt MoveCount)
{
    m_Exponent-=MoveCount;
    TMaxInt i=MoveCount/em10Power;
    TMaxInt iMoveCount= MoveCount%em10Power;
    
   //完成大的移动
 int m_CalcSize=m_Digits.size();
    if (i>=m_CalcSize)
    {
        Clear();
        return;
    }
    else if (i>0) //左移i
    {
        for (int j=0;j<m_CalcSize-i;++j)
            m_Digits[j]=m_Digits[int(j+i)];
        for (int k=int(m_CalcSize-i);k<m_CalcSize;++k)
            m_Digits[k]=0;
    }
    
    if (iMoveCount>0) //完成小的移动
    {
        int iMul=1;
        for (int k=0;k<iMoveCount;++k) iMul*=10;

  for (int i=m_CalcSize-1;i>=0;--i)
  {
   m_Digits[i]*=iMul;
  }
  for (int j=m_CalcSize-1;j>=1;--j)
  {
   if (m_Digits[j]>emBase)
   {
    m_Digits[j-1]+=m_Digits[j]/emBase;
    m_Digits[j]%=emBase;
   }
  }
    }
}

void TLargeFloat::MoveRight10Power(TMaxInt MoveCount)
{
    m_Exponent+=MoveCount;
    TMaxInt i=MoveCount/em10Power;
    TMaxInt iMoveCount= MoveCount%em10Power;
    
   //完成大的移动
 int m_CalcSize=m_Digits.size();
    if (i>=m_CalcSize)
    {
        Clear();
        return;
    }
    else if (i>0) //右移i
    {
        for (int j=m_CalcSize-1;j>=i;--j)
            m_Digits[j]=m_Digits[int(j-i)];
        for (int k=0;k<i;++k)
            m_Digits[k]=0;
    }
    if (iMoveCount>0) //完成小的移动
    {
        Int32bit iDiv=1;
        for (int j=0;j<iMoveCount;++j) iDiv*=10;
  
  for ( i=0;i<m_CalcSize-1;++i)
  {
   m_Digits[int(i)+1]+=(m_Digits[int(i)]%iDiv)*emBase;
   m_Digits[int(i)]/=iDiv;
  }
  if (m_CalcSize>=1)
   m_Digits[m_CalcSize-1]/=iDiv;
    }
}

void TLargeFloat::Canonicity()
{
    //规格化 小数部分
 //规格化的数满足的条件：小数部分的值 小于1.0，大于等于0.1;

    if (m_Sign==0)
    {
        Clear();
        return;
    }

    int size=m_Digits.size();
    if (m_Digits[0]>=emBase)//是否需要右移1
    {
        m_Exponent+=em10Power;
        for (int j=size-1;j>=2;--j)
            m_Digits[j]=m_Digits[j-1];
        if (size>=2) m_Digits[1]=m_Digits[0] % emBase;
        m_Digits[0]/=emBase;
    }
 

    //考虑左移
    int i=0;
    int iMoveCount=0;
    for (;i<size;++i)//寻找第一个不为零的位i
    {
        if (m_Digits[i]>0)
  {
   Int32bit iBase=emBase;//iMoveCount 需要移动的十进制位数
   for (;iMoveCount<em10Power;++iMoveCount)
   {
    iBase/=10;
    if (m_Digits[i]>=iBase) break;
   }
   break;
  }
    }
    MoveLeft10Power(i*em10Power+iMoveCount);

}

std::string  TLargeFloat::DigitToString(Int32bit iDigit)
{
    char buffer[40];
 for (int j=0;j<40;++j) buffer[j]=0;
    _itoa(iDigit, buffer,10);
    std::string result(buffer);
 assert(em10Power>=result.size());
    result=std::string(em10Power-result.size(),'0')+result;
    return result;
}

std::string  TLargeFloat::AsString() const//转换为字符串
{
 //没有优化

    if (m_Sign==0)
        return "0";
    else 
    {
        std::string result;
        if  (m_Sign<0)
            result="-0.";
        else
            result="0.";
        int size=m_Digits.size();
        for (int i=0;i<size;++i)
        {
            result+=DigitToString(m_Digits[i]);
        }

        if (m_Exponent.AsInt()!=0)
        {
            result+='e';
            if (m_Exponent.AsInt()<0) result+='-';
            char buffer[40];
   //_i64toa(Private_::abs(m_Exponent.AsInt()), buffer,10);//__int64
   _itoa(Private_::abs(m_Exponent.AsInt()), buffer,10);
            result+=buffer;
        }
        return result;
    }
}


void TLargeFloat::Chs()//求负
{
    m_Sign*=(-1);
}

void TLargeFloat::Abs()//绝对值
{
    if (m_Sign!=0) m_Sign=1;
}


void TLargeFloat::swap(SelfType& Value)
{
    m_Digits.swap(Value.m_Digits);
 std::swap(m_Sign,Value.m_Sign);
    std::swap(m_Exponent,Value.m_Exponent);
}

//////////////
void TLargeFloat::toEqExponent(SelfType& x,SelfType& y)//对齐小数点
{
 ExpInt MoveCount=x.m_Exponent;
 MoveCount-=y.m_Exponent;
 if (MoveCount.AsInt()<0) MoveCount=ExpInt(-MoveCount.AsInt());
    if (x.m_Exponent.AsInt()>y.m_Exponent.AsInt())
    {
  y.MoveRight10Power(MoveCount.AsInt());
        y.m_Exponent=x.m_Exponent;
    }
    else if(x.m_Exponent.AsInt()<y.m_Exponent.AsInt())
    {
  x.MoveRight10Power(MoveCount.AsInt());
        x.m_Exponent=y.m_Exponent;
    }
}

void TLargeFloat::Abs_Add(const SelfType& Value)//绝对值加法
{
    SelfType Right(Value);
    Abs();
    Right.Abs();

    toEqExponent(Right,*this);//小数点对齐
    SetSameSizeMax(Right,*this);

 int m_CalcSize=m_Digits.size();
 
    m_Digits[0]+=Right.m_Digits[0];
    for (int j=m_CalcSize-1;j>=1;--j)
    {
        m_Digits[j]+=Right.m_Digits[j];
        if (m_Digits[j]>=emBase)//进位
        {
            m_Digits[j-1]+=1;
            m_Digits[j]-=emBase;
        }
    }

 if (m_Digits[0]>=emBase)
  Canonicity();

}

void TLargeFloat::Abs_Sub_Abs(const SelfType& Value)//绝对值减法
{

    SelfType Right(Value);
    Abs();
    Right.Abs();


    toEqExponent(Right,*this);//小数点对齐
    SetSameSizeMax(Right,*this);

    int comResult=Compare(Right);

    if (comResult==0)
    {
        Clear();
        return;
    }
    
 if (comResult<0)
        swap(Right);

 int m_CalcSize=m_Digits.size();
 for (int j=m_CalcSize-1;j>=0;--j)
    {
        m_Digits[j]-=Right.m_Digits[j];
  
  if (m_Digits[j]<0)//借位
  {
   m_Digits[j-1]-=1;//j-1不可能超界
   m_Digits[j]+=emBase;
  }
    
 }

 if (m_Digits[0]*10<emBase)
  Canonicity();
}

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

TLargeFloat::SelfType& TLargeFloat::operator += (const SelfType& Value)
{
 if (Value.m_Sign==0)
  return (*this);
 else if (m_Sign==0)
  return (*this)=Value;
 else if (m_Sign==Value.m_Sign)
    {
        Int32bit oldSign=m_Sign;
        Abs_Add(Value);
        m_Sign=oldSign;
    }
    else
    {
        SelfType x(Value);
        x.Chs();
        *this-=(x);
    }
    return *this;
}

TLargeFloat::SelfType& TLargeFloat::operator -= (const SelfType& Value)
{
 if (Value.m_Sign==0)
  return (*this);
 else if (m_Sign==0)
 {
  (*this)=Value;
  this->Chs();
  return (*this);
 }
    else if (m_Sign==Value.m_Sign)
    {
        int comResult=Compare(Value);    
        if (comResult==0)
        {
            Clear();
        }
        else
        {
            Abs_Sub_Abs(Value);
            m_Sign=comResult;
        } 
    }
    else
    {
        SelfType x(Value);
        x.Chs();
        *this+=(x);
    }
    return *this;
}


int  TLargeFloat::Compare(const SelfType& Value) const//*this>Value 返回1，小于返回-1，相等返回0
{
    //
    if (m_Sign>Value.m_Sign)
        return 1;
    else if(m_Sign<Value.m_Sign)
        return -1;
    else if(m_Sign==0)//m_Sign==Value.m_Sign
        return 0;

    //m_Sign==Value.m_Sign
    if (m_Exponent.AsInt()>Value.m_Exponent.AsInt())
    {
        if (m_Sign<0) 
            return -1;
        else
            return 1;
    }
    else if (m_Exponent.AsInt()<Value.m_Exponent.AsInt())
    {
        if (m_Sign<0) 
            return 1;
        else
            return -1;
    }
    else//(m_Exponent==Value.m_Exponent)
    {
        int sizeS=m_Digits.size();
        int sizeV=Value.m_Digits.size();
        int size=Private_::min(sizeS,sizeV);
        int result=0;
        for (int i=0;i<size;++i)
        {
            if (m_Digits[i]>Value.m_Digits[i])
            {
                result=1;
                break;
            }
            else if (m_Digits[i]<Value.m_Digits[i])
            {
                result=-1;
                break;
            }
        }
        if (result==0)
        {   //继续比较 处理尾部
            if (sizeS>sizeV)
            {
                for (int i=sizeV;i<sizeS;++i) 
                {
                    if (m_Digits[i]>0)
                    {
                        result=1;
                        break;
                    }
                }
            }
            else if (sizeS<sizeV)
            {
                for (int i=sizeS;i<sizeV;++i) 
                {
                    if (Value.m_Digits[i]>0)
                    {
                        result=-1;
                        break;
                    }
                }
            }
        }
        return result*m_Sign;
        //
    }

    
}

void TLargeFloat::SetSameSizeMax(SelfType& x,SelfType& y)//调整有效位数
{
    int sizeX=x.m_Digits.size();
    int sizeY=y.m_Digits.size();
    if (sizeX<sizeY)
    {
        x.m_Digits.resize(sizeY,0);
    }
    else if (sizeX>sizeY)
    {
        y.m_Digits.resize(sizeX,0);
    }
}

TLargeFloat::SelfType& TLargeFloat::operator *= (const SelfType& Value)
{
    if ((Value.m_Sign==0)||(m_Sign==0)) 
    {
        Clear();
        return *this;
    }

 //先考虑符号和指数
    int sign=Value.m_Sign*m_Sign;
    ExpInt oldExponent=Value.m_Exponent;
 oldExponent+=m_Exponent;

    SelfType x(Value);
    SelfType y(*this); 
    x.Abs();
    y.Abs();
    x.m_Exponent=0;
    y.m_Exponent=0;

    SetSameSizeMax(x,y);
 unsigned int mulSize=x.m_Digits.size();

 TArray tempArray(mulSize*2);
 ArrayMUL(&(x.m_Digits[0]),&(y.m_Digits[0]),&(tempArray[0]),mulSize);
 tempArray.resize(x.m_Digits.size());
 this->m_Digits.swap(tempArray);

    m_Sign=sign;
    m_Exponent=-em10Power;
    m_Exponent+=(oldExponent);
    Canonicity();
    return *this;
}

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
 // 证毕.

    ExpInt oldExponent=m_Exponent;
    m_Exponent=0;

    SelfType& a=*this;

    SelfType x(a);
 long double fx=x.AsFloat();
 if (fx==0) throw TException("ERROR:TLargeFloat::Rev()");

    x=1.0/fx;//初始迭代值
    
    SelfType temp(0.0,TLargeFloat::TDigits(x.GetDigitsLength()));
    int size=x.m_Digits.size();
    int n=6;//n用来用来追踪计算出的有效精度,并防止死循环

    while (true)
    {
  temp=(2-a*x)*x;
  if (temp.Compare(x)==0) break;
        x=temp;

  if (n>size) break;
  n*=2;//二阶收敛
    }

    x.m_Exponent-=oldExponent;
    swap(x);
}