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

//TLargeFloat.cpp
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


main函數，Demo程序
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// LargeFloat.cpp : Defines the entry point for the console application.
//

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

int main(int argc, char* argv[])
{
 
 try
 {
  LargeFloat_UnitTest();
  
  TLargeFloat x(0,TLargeFloat::TDigits(100));
  x=1;
  Debug_toCout("1/7:=",x/7);
  
  x=double(123456787654321);// __int64(123456787654321)
  Debug_toCout("123456787654321*123456787654321 := ",x*x);

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

  //test PI
  std::cout<<std::endl;
  x=GetBorweinPI();
  Debug_toCout("PI := ",x);

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

}
/////////////////////