tlargefloat.cpp

利用数值分析方法

            }
        }
        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);
}


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阶收敛公式
 // 证明略
    ExpInt sqrExponent=m_Exponent.AsInt()/2;
    m_Exponent-=sqrExponent;
    m_Exponent-=sqrExponent;

    SelfType& a=*this;

    SelfType x(a);
  long double fx=x.AsFloat();
 if (fx<=0) throw TException("ERROR:TLargeFloat::RevSqrt()");
 x=1.0/sqrt(fx);//初始迭代值
    
    SelfType temp(0.0,TLargeFloat::TDigits(x.GetDigitsLength()));
    int size=x.m_Digits.size();
    int n=6;//n用来用来追踪计算出的有效精度,并防止死循环
    while (true)
    {
        temp=(3-a*x*x)*x/2;
  if (temp.Compare(x)==0) break;
        x=temp;

  if (n>size) break;
  n*=2;//2阶收敛公式
    }
    x.m_Exponent-=sqrExponent;
    swap(x);
    
}

void TLargeFloat::Sqrt() //求x^0.5;
{
    SelfType x(*this);
    x.RevSqrt();
    (*this)*=x;
}

void TLargeFloat::ArrayMUL(const Int32bit* x,const Int32bit* y,Int32bit* result, unsigned int MulSize)
{

 //没有优化 

 for (int k=0;k<MulSize*2;++k) result[k]=0;

    for (int i=0;i<int(MulSize);++i)
    {
        if (x[i]>0)
        {
            for (int j=0;j<int(MulSize);++j)
            {
                int index=i+j;
                result[index]+=(x[i]*y[j]);
            }
            for (int k=MulSize*2-1;k>=1;--k)
            {
                if (result[k]>=emBase)
                {
                    result[k-1]+=result[k]/emBase;
                    result[k]=result[k] % emBase;
                }
            }
        }
    }
}

void TLargeFloat::MulInt(TMaxInt iValue)//乘以一个整数
{
 if (iValue==0) 
 {
  Clear();
  return;
 }
 else if (iValue<0)
 {
  iValue=-iValue;
  this->Chs();
  //continue
 }

 Int32bit iValueLow=Int32bit(iValue % emBase);
 TMaxInt iValueHigh=iValue / emBase;

 TLargeFloat temp(*this);

 int size=m_Digits.size();
 for (int i=size-1;i>=0;--i)
 {
  m_Digits[i]*=iValueLow;
 }
 for (int j=size-1;j>=1;--j)
 {
  if (m_Digits[j]>emBase)
  {
   m_Digits[j-1]+=m_Digits[j]/emBase;
   m_Digits[j]%=emBase;
  }
 }
 Canonicity();

 if (iValueHigh!=0)
 {
  temp.m_Exponent+=em10Power;
  temp.MulInt(iValueHigh);//递归
  (*this)+=temp;
 }
}

void TLargeFloat::DivInt(TMaxInt iValue)//除以一个整数
{
 if (iValue==0) 
 {
  throw TException("ERROR:TLargeFloat::DivInt()");
 }
 else if (iValue<0)
 {
  iValue=-iValue;
  this->Chs();
  //continue
 }

 if (iValue<=emBase)
 {
  Int32bit iDiv=Int32bit(iValue);
  int size=m_Digits.size();
  for (int i=0;i<size-1;++i)
  {
   m_Digits[i+1]+=(m_Digits[i]%iDiv)*emBase;
   m_Digits[i]/=iDiv;
  }
  if (size>=1)
   m_Digits[size-1]/=iDiv;
  Canonicity();
 }
 else
 {
  int size=m_Digits.size();
  TMaxInt HighData=0;
  for (int i=0;i<size-1;++i)
  {
   TMaxInt tempData=(HighData*emBase+m_Digits[i])%iValue;
   m_Digits[i]=Int32bit((HighData*emBase+m_Digits[i])/iValue);
   HighData=tempData;
  }
  if (size>=1)
   m_Digits[size-1]=Int32bit((HighData*emBase+m_Digits[size-1])/iValue);
  Canonicity();
 }

}


const TLargeFloat operator + (const TLargeFloat& x,const TLargeFloat& y)
{
 TLargeFloat temp(x);
 return temp+=y;
}

 const TLargeFloat operator - (const TLargeFloat& x,const TLargeFloat& y)
{
 TLargeFloat temp(x);
 return temp-=y;
}
 const TLargeFloat operator * (const TLargeFloat& x,const TLargeFloat& y)
{
 TLargeFloat temp(x);
 return temp*=y;
}
 const TLargeFloat operator / (const TLargeFloat& x,const TLargeFloat& y)
{
 TLargeFloat temp(x);
 return temp/=y;
}


 const TLargeFloat::SelfType TLargeFloat::operator -  () const//求负
{
 SelfType temp(*this);
 temp.Chs();
 return temp;
}

 const TLargeFloat::SelfType& TLargeFloat::operator +  () const//求正
{
 return *this;
}

TLargeFloat::SelfType& TLargeFloat::operator *= (long double  fValue)
{
 if (!FloatIsInteger(fValue))
 {
  TLargeFloat temp(fValue,TCharacter());
  return (*this)*=temp;
 }
 else
 {
  TMaxInt iValue=TMaxInt(floor(fValue));
  MulInt(iValue);
  return *this;
 }
}

TLargeFloat::SelfType& TLargeFloat::operator /= (long double  fValue)
{
 if (!FloatIsInteger(fValue))
 {
  TLargeFloat temp(fValue,TCharacter());
  return (*this)/=temp;
 }
 else
 {
  TMaxInt iValue=TMaxInt(floor(fValue));
  DivInt(iValue);
  return *this;
 }
}

 TLargeFloat::SelfType& TLargeFloat::operator += (long double  fValue)
{
 TLargeFloat temp(fValue,TCharacter());
 return (*this)+=temp;
}

 TLargeFloat::SelfType& TLargeFloat::operator -= (long double  fValue)
{
 TLargeFloat temp(fValue,TCharacter());
 return (*this)-=temp;
}

 const TLargeFloat operator + (const TLargeFloat& x,long double y)
{
 TLargeFloat temp(x);
 return temp+=y;
}
 const TLargeFloat operator - (const TLargeFloat& x,long double y)
{
 return x+(-y);
}
 const TLargeFloat operator * (const TLargeFloat& x,long double y)
{
 TLargeFloat temp(x);
 return temp*=y;
}
 const TLargeFloat operator / (const TLargeFloat& x,long double y)
{
 TLargeFloat temp(x);
 return temp/=y;
}
 const TLargeFloat operator + (long double x,const TLargeFloat& y)
{
 return y+x;
}
 const TLargeFloat operator - (long double x,const TLargeFloat& y)
{
 return (-y+x);
}
 const TLargeFloat operator * (long double x,const TLargeFloat& y)
{
 return y*x;
}
 const TLargeFloat operator / (long double x,const TLargeFloat& y)
{
 TLargeFloat temp(y);
 temp.Rev();
 temp*=x;
 return temp;
}

 bool operator ==(const TLargeFloat& x,const TLargeFloat& y)
{
 return (x.Compare(y)==0);
}

 bool operator < (const TLargeFloat& x,const TLargeFloat& y)
{
 return (x.Compare(y)<0);
}

 bool operator ==(const TLargeFloat& x,long double y)
{
 TLargeFloat tempy(y,TLargeFloat::TCharacter());
 return (x==tempy);
}
 bool operator < (const TLargeFloat& x,long double y)
{
 TLargeFloat tempy(y,TLargeFloat::TCharacter());
 return (x<tempy);
}


/////////////////////
//UnitTest
//
namespace Private_UnitTest
{

#define _TestCheck(Value) \
{\
 if (!(Value))\
 {\
  throw TLargeFloatException("UnitTest not pass when Check:\"" #Value "\""); \
 }\
}

#define _TestCheck_Exception(Value,LabLine)/*检测异常*/ \
{\
 int __TestCheck_Exception_ERROR =0; \
 try\
 {\
  Value; \
 }\
    catch(TLargeFloatException&)\
 {\
  goto ok_lab_##LabLine;\
 } \
 throw TLargeFloatException("_TestCheck_Exception not pass when Check:\"" #Value "\""); \
    ok_lab_##LabLine: ;\
}

void LargeFloat_BaseTest()
{
 TLargeFloat x(0,TLargeFloat::TDigits(100));
 _TestCheck(x==0);
 
 x=1.0e20;
 x=-x;

 TLargeFloat y=x;
 _TestCheck(y==-1.0e20);
 _TestCheck(x==y);

 x=1.0e-4;
 y=1.0e-11;
 TLargeFloat z=x;
 z+=y;
 _TestCheck(z==(x+y));

 z=y-x;
 _TestCheck(z+x==y);

 const long double ttmpld=0.000123454365465454765474;
 x=ttmpld;
 _TestCheck(x==ttmpld);
 _TestCheck(Private_::abs(x.AsFloat()-ttmpld)<1e-10);


 //const __int64 ttmpi64=-1234567891235555;
 const double ttmpi64=-1234567891235555;
 x=ttmpi64;
 _TestCheck(x==ttmpi64);
 _TestCheck(abs(x/ttmpi64-1)<0.00001);
 _TestCheck(abs(x*x-x*ttmpi64)<0.00001);
 _TestCheck(abs(x.AsFloat()-ttmpi64)<0.00001);

 x=-x;
 _TestCheck(x==-ttmpi64);
 x=0.0;
 _TestCheck(0==x);

 const int tmpi=-123;
 x=tmpi;
 _TestCheck(-tmpi==abs(x));
 _TestCheck(tmpi*tmpi==sqr(x));

}

void LargeFloat_TestException()
{
 TLargeFloat  x(-2,TLargeFloat::TDigits(1000));
 TLargeFloat  y(-2,TLargeFloat::TDigits(1000));

 _TestCheck_Exception(sqrt(x),0);
 _TestCheck_Exception(x/0,1);
 _TestCheck_Exception(x/(y=0),2);
 /*
 x=1.1234e300;
 x=x*x;
 _TestCheck_Exception(x.AsFloat(); ,3);

 x=1.1234e150;
 _TestCheck_Exception( for(;;) {x=x*x;} ,4);
 x=1.1234e-200;
 _TestCheck_Exception( for(;;) {x=x*x;} ,5);
  */

 _TestCheck_Exception(TLargeFloat(545,TLargeFloat::TDigits(0)); , 6);

 x=0;
 _TestCheck_Exception( revsqrt(x); ,7);

}

void LargeFloat_TestCompare()
{
 TLargeFloat  x(-0.3456346e-5,TLargeFloat::TDigits(1000));
 TLargeFloat  y(3.5453,TLargeFloat::TDigits(1000));

 _TestCheck(x!=y);
 _TestCheck(x<y);
 _TestCheck(y>x);

 y=x;
 _TestCheck(x==y);
 _TestCheck(!(x<y));
 _TestCheck(!(y>x));

 _TestCheck(x!=0);
 _TestCheck(x<0);
 _TestCheck(0>x);
 x=123;
 _TestCheck(x==123);
 _TestCheck(!(x<123));
 _TestCheck(!(123>x));
}

void LargeFloat_TestAbs_Add()
{
 //todo:

}

//todo: 其他测试

}//end namespace Private_UnitTest


void LargeFloat_UnitTest()//正确性测试
{
 using namespace Private_UnitTest;

 LargeFloat_BaseTest();
 LargeFloat_TestException();
 LargeFloat_TestCompare();
 LargeFloat_TestAbs_Add();
}

////

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