自己动手打造“超高精度浮点数类”
自己动手打造“超高精度浮点数类”
很多人可能都想自己写一个能够执行任意精度计算的浮点数；:D我写的第一个程序就是用qbasic计算自然数e到100万位(后来计算PI); 这里有一个C++类的实现TLargeFloat :
（共有3个文件；在vc6和dev-c++中编译通过；里面有一个计算PI的Borwein四次迭代式）
类的声明文件:TLargeFloat.h
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// TLargeFloat.h: interface for the TLargeFloat class.
// 超高精度浮点数类TLargeFloat
// 2004.03.28 by HouSisong@263.net
//////////////////////////////////////////////////////////////////////
// Update 2004.03.31 by HouSisong
// Update 2004.04.05 by HouSisong
// Update 2004.04.13 by HouSisong 添加异常触发能力,TLargeFloat捕获所有非法运算抛出TLargeFloatException异常
#ifndef _TLARGE_FLOAT_H__INCLUDED_
#define _TLARGE_FLOAT_H__INCLUDED_
#include <vector>
#include <sstream>
#include <string>
#include <Exception>
#include <limits>
#include <algorithm>
class TLargeFloat;//超高精度浮点数类TLargeFloat
class TLargeFloatException;//超高精度浮点数异常类
//改进方向：
// 1.强力优化ArrayMUL数组乘运算(当前算法复杂度为n*n)，
// 可以使用二分算法来降低运算量，并使用复杂度为n*log(n)的快速复利叶变换或数论变换）
// 2.增加运算精度动态控制能力，有利于优化，减少乘法量；
// 3.添加新的基本运算函数，如：指数运算power、对数运算log、三角函数sin,cos,tan等
// 4.可以考虑：内部使用2的次方的底数;这样的话，输出函数就会麻烦一些了
//////注意：如果浮点数与TLargeFloat进行混合运算；
// 可能会产生误差(有效位数会受到浮点数影响)；
// 整数 或 为可表示整数的浮点数 参与运算不会产生误差；
//超高精度浮点数异常类
class TLargeFloatException :public std::exception
{
private:
std::string m_ErrorMsg;
public:
TLargeFloatException() {};
TLargeFloatException(const char * Error_Msg) :m_ErrorMsg(Error_Msg){ }
virtual const char* what() const throw() { return m_ErrorMsg.c_str();}
virtual ~TLargeFloatException() throw() {}
};
//TCatchIntError只是对整数类型TInt进行的包装
//设计TCatchIntError是为了当整数运算超出值域的时候，抛出异常
//超高精度浮点数类的指数运算时使用
template <typename TInt,typename TException,TInt MinValue,TInt MaxValue>
//<要包装的整数类型，超界时抛出的异常类型，TInt最小值，TInt最大值>
class TCatchIntError
{
private:
typedef TCatchIntError<TInt,TException,MinValue,MaxValue> SelfType;
TInt m_Int;
SelfType& inc(TInt uValue)
{
if (MaxValue-uValue<m_Int)
throw TException("ERROR:TCatchIntError::inc(); ");
m_Int+=uValue;
return (*this);
}
SelfType& dec(TInt uValue)
{
if (MinValue+uValue>m_Int)
throw TException("ERROR:TCatchIntError::dec()");
m_Int-=uValue;
return (*this);
}
public:
TCatchIntError() :m_Int(0){}
TCatchIntError(TInt Value) :m_Int(Value){}
TCatchIntError(const SelfType& Value) :m_Int(Value.m_Int){}
TInt AsInt()const { return m_Int; }
SelfType& operator +=(TInt Value) //throw(TLargeFloatException)
{ if (Value<0) return dec(-Value);
else return inc(Value); }
SelfType& operator -=(TInt Value) //throw(TLargeFloatException)
{ if (Value<0) return inc(-Value);
else return dec(Value); }
SelfType& operator +=(const SelfType& Value) { return (*this)+=(Value.m_Int); }//throw(TLargeFloatException)
SelfType& operator -=(const SelfType& Value) { return (*this)-=(Value.m_Int); }//throw(TLargeFloatException)
};
////填写编译器支持的较大的整数类型
//__int64 Int64_Min() { return std::numeric_limits<__int64>::min(); }//返回0, :(
//__int64 Int64_Max() { return std::numeric_limits<__int64>::max(); }//返回0, :(
//const __int64 Int64_Min = - __int64(9223372036854775808);//注意负号
//const __int64 Int64_Max = __int64(9223372036854775807);
const long int Int64_Min = -2147483648;//注意负号
const long int Int64_Max = 2147483647;
namespace Private_
{
template<typename T>
inline const T& min(const T& x,const T& y)//求最小值
{
if (x>y)
return y;
else
return x;
}
template<typename T>
inline const T& max(const T& x,const T& y)//求最大值
{
if (x>y)
return x;
else
return y;
}
template<typename T>
inline const T abs(const T& x)//求绝对值
{
if (x<0)
return -x;
else
return x;
}
template<typename T>
inline void swap(T& x,T& y) //交换两个变量的值
{
T temp=x;
x=y;
y=temp;
}
}//end namespace
//超高精度浮点数类
class TLargeFloat
{
private:
enum {
emLongDoubleDigits=std::numeric_limits<long double>::digits10,//long double的10进制有效精度
emLongDoubleMaxExponent=std::numeric_limits<long double>::max_exponent10,//long double的最大10进制指数
emLongDoubleMinExponent=std::numeric_limits<long double>::min_exponent10 };//long double的最小10进制指数
typedef TLargeFloat SelfType;
typedef TLargeFloatException TException;
typedef long int Int32bit;//32bit位的整数类型,超过也可以
//typedef __int64 TMaxInt; //填写编译器支持的较大的整数类型
typedef long int TMaxInt; //填写编译器支持的较大的整数类型
typedef TCatchIntError<TMaxInt,TException,Int64_Min,Int64_Max> ExpInt;//注意：后面的两个值为TMaxInt的最小值和最大值
typedef std::vector<Int32bit> TArray;//小数位使用的数组类型
enum { em10Power=4, emBase=10000};//数组为10000进制，数组的一个元素表示一位，对应4个十进制位
Int32bit m_Sign; //符号位 正:1, 负:-1, 零: 0
ExpInt m_Exponent; //保存10为底的指数
TArray m_Digits; //小数部分 排列顺序是TArray[0]为第一个小数位,依此类推；取值范围0--999
void Abs_Add(const SelfType& Value);//绝对值加 x:=|x|+|y|;
void Abs_Sub_Abs(const SelfType& Value);//绝对值减的绝对值x:=| |x|-|y| |;
void MoveLeft10Power(TMaxInt MoveCount);//十进制指数移动 值不变指数增大MoveCount
void MoveRight10Power(TMaxInt MoveCount);//十进制指数移动 值不变指数减小MoveCount
void MulInt(TMaxInt iValue);//乘以一个整数；
void DivInt(TMaxInt iValue);//除以一个整数；
void Clear();//清零
void Chs();//求负
int Compare(const SelfType& Value) const;//比较两个数；(*this)>Value 返回1，小于返回-1，相等返回0
void Canonicity();//规格化 转化值到合法格式
static std::string DigitToString(Int32bit iDigit);//将数组的一个元素转换为字符串表示
static void toEqExponent(SelfType& x,SelfType& y);//值不变，x,y的小数点对齐
static void SetSameSizeMax(SelfType& x,SelfType& y);//使两个高精度数的有效位数相同，位数小的进行提升
static bool FloatIsInteger(long double fValue);//判断浮点数是否为可表示整数
static unsigned int DigitsSize(unsigned int uiDigitsLength);//
//数组乘 （卷积result[i+j]=x[i]*y[j];) //ArrayMUL 是需要优化的首要目标
static void ArrayMUL(const Int32bit* x,const Int32bit* y,Int32bit* result,unsigned int MulSize);
class TCharacter{};
TLargeFloat(long double DefultFloatValue,const TCharacter&);//内部使用，浮点数转化为 TLargeFloat
void Abs();//绝对值
void Rev();//求倒数1/x
void RevSqrt();//求1/x^0.5;
void Sqrt();//求x^0.5;
public:
class TDigits//TDigits用来设置TLargeFloat的精度;//增加这个类是为了避免TLargeFloat的构造函数的可能误用
{
private:
unsigned int m_eDigits;
public:
explicit TDigits(unsigned int uiDigitsLength) :m_eDigits(uiDigitsLength){}
unsigned int GetDigits()const { return m_eDigits; }
};
TLargeFloat(const SelfType& Value);
TLargeFloat(long double DefultValue,const TDigits& DigitsLength);//TDigits (十进制的)有效位数
virtual ~TLargeFloat(){}
void swap(SelfType& Value);//交换值
unsigned int GetDigitsLength() const;//返回当前的10进制有效位数
void SetDigitsLength(unsigned int uiDigitsLength);//重新设置10进制有效位数
void SetDigitsLength(const TDigits& DigitsLength) { SetDigitsLength(DigitsLength.GetDigits()); }
long double AsFloat() const;//转化为浮点数
std::string AsString() const;//转换为字符串
const SelfType operator - () const;//求负 //注意：不能使用SelfType&
const SelfType& operator + () const;//求正 //注意：可以使用SelfType&,因为值不变
SelfType& operator = (long double fValue); //注意：转换可能存在小的误差
SelfType& operator = (const SelfType& Value); //编译器默认的也行
SelfType& operator *= (long double fValue);
SelfType& operator /= (long double fValue);
SelfType& operator += (long double fValue);
SelfType& operator -= (long double fValue);
SelfType& operator += (const SelfType& Value);
SelfType& operator -= (const SelfType& Value);
SelfType& operator *= (const SelfType& Value);
SelfType& operator /= (const SelfType& Value);
friend const TLargeFloat operator + (const TLargeFloat& x,const TLargeFloat& y);
friend const TLargeFloat operator - (const TLargeFloat& x,const TLargeFloat& y);
friend const TLargeFloat operator * (const TLargeFloat& x,const TLargeFloat& y);
friend const TLargeFloat operator / (const TLargeFloat& x,const TLargeFloat& y);
friend const TLargeFloat operator + (const TLargeFloat& x,long double y);
friend const TLargeFloat operator - (const TLargeFloat& x,long double y);
friend const TLargeFloat operator * (const TLargeFloat& x,long double y);
friend const TLargeFloat operator / (const TLargeFloat& x,long double y);
friend const TLargeFloat operator + (long double x,const TLargeFloat& y);
friend const TLargeFloat operator - (long double x,const TLargeFloat& y);
friend const TLargeFloat operator * (long double x,const TLargeFloat& y);
friend const TLargeFloat operator / (long double x,const TLargeFloat& y);
friend bool operator ==(const TLargeFloat& x,const TLargeFloat& y);
friend bool operator < (const TLargeFloat& x,const TLargeFloat& y);
friend bool operator ==(const TLargeFloat& x,long double y);
friend bool operator < (const TLargeFloat& x,long double y);
friend bool operator ==(long double x,const TLargeFloat& y) { return (y==x); }
friend bool operator < (long double x,const TLargeFloat& y) { return (y>x); }
friend bool operator !=(const TLargeFloat& x,const TLargeFloat& y) { return !(x==y); }
friend bool operator > (const TLargeFloat& x,const TLargeFloat& y) { return (y<x); }
friend bool operator >=(const TLargeFloat& x,const TLargeFloat& y) { return !(x<y); }
friend bool operator <=(const TLargeFloat& x,const TLargeFloat& y) { return !(x>y); }
friend bool operator !=(const TLargeFloat& x,long double y) { return !(x==y); }
friend bool operator > (const TLargeFloat& x,long double y) { return (y<x); }
friend bool operator >=(const TLargeFloat& x,long double y) { return !(x<y); }
friend bool operator <=(const TLargeFloat& x,long double y) { return !(x>y); }
friend bool operator !=(long double x,const TLargeFloat& y) { return !(x==y); }
friend bool operator > (long double x,const TLargeFloat& y) { return (y<x); }
friend bool operator >=(long double x,const TLargeFloat& y) { return !(x<y); }
friend bool operator <=(long double x,const TLargeFloat& y) { return !(x>y); }
friend std::ostream& operator << (std::ostream& cout, const TLargeFloat& Value) { return cout<<Value.AsString(); }
friend void swap(TLargeFloat& x,TLargeFloat& y) { x.swap(y); }
friend const TLargeFloat abs(const TLargeFloat& x) { TLargeFloat result(x); result.Abs(); return result; }//绝对值,|x|
friend const TLargeFloat sqrt(const TLargeFloat& x) { TLargeFloat result(x); result.Sqrt(); return result;} //开方,x^0.5
friend const TLargeFloat revsqrt(const TLargeFloat& x) { TLargeFloat result(x); result.RevSqrt(); return result; }//求1/x^0.5;
friend const TLargeFloat sqr(const TLargeFloat& x) { return x*x; };//平方,x^2
};
void LargeFloat_UnitTest();//正确性测试
//下面的代码是用来测试的
//例子：
//计算圆周率PI
//经过测试,计算中97%的时间都在运行ArrayMUL函数
TLargeFloat GetBorweinPI();//使用Borwein四次迭代式
void Debug_toCout(const std::string& strx,const TLargeFloat& x);//调试输出
#endif // _TLARGE_FLOAT_H__INCLUDED_
// TLargeFloat.h
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
类的实现文件:TLargeFloat.cpp
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// TLargeFloat.cpp: implementation of the TLargeFloat class.
// 超高精度浮点数类TLargeFloat
//////////////////////////////////////////////////////////////////////
//#include "stdafx.h"//
#include "TLargeFloat.h"
#include "assert.h"
#include <math.h>
#include <iostream>
//计算圆周率PI
TLargeFloat GetBorweinPI()
{
/*
Borwein四次迭代式：
初值：
a0=6-4*2^0.5; y0=2^0.5-1;
重复计算：
y(n+1)=[1-(1-y^4)^0.25]/[1+(1-y^4)^0.25];
y=y(n+1);
a(n+1)=a*(1+y)^4-2^(2*n+3)*y*(1+y+y*y);
最后计算：
PI=1/a;
*/
TLargeFloat::TDigits sCount(1000);//计算采用1000位精度 计算出的PI后面35位无效
TLargeFloat a(0.0,sCount),y(0.0,sCount);
TLargeFloat PI(0.0,sCount);
TLargeFloat temp(0.0,sCount),temp2(0.0,sCount);
TLargeFloat pow(0.0,sCount);
pow=(2*2*2);//2^(2*n+3); (n=0);
//a0=6-4*2^0.5;
temp=2;
temp=sqrt(temp);
a=6-4*temp;
//y0=2^0.5-1;
y=temp-1;
//1
int m=8;//
int n=0;
while (true)
{
//y(n+1)=[1-(1-y^4)^0.25]/[1+(1-y^4)^0.25];
temp=1-sqr(sqr(y));
temp=revsqrt(revsqrt(temp));//等价于 temp=sqrt(sqrt(temp));
y=(1-temp)/(1+temp);
//a(n+1)=a*(1+y)^4-2^(2*n+3)*y*(1+y+y*y);
temp=sqr(sqr(1+y));
a=a*temp-pow*y*(1+y+y*y);
pow*=4;
++n;
if (m>int(sCount.GetDigits())) break;
m*=4;//四阶收敛
}
//PI=1/a;
PI=1/a;
return PI;
}
/////////
bool TLargeFloat::FloatIsInteger(long double fValue)//浮点数是否为可表示整数
{
return (TMaxInt(floor(fValue))==fValue);
}
unsigned int TLargeFloat::DigitsSize(unsigned int uiDigitsLength)
{
if (!(uiDigitsLength>=1))
{
throw TException("ERROR:TLargeFloat::DigitsSize()");
}
return (uiDigitsLength+(em10Power-1))/em10Power;//计算出需要的数组大小
}
TLargeFloat::TLargeFloat(long double DefultFloatValue,const TDigits& DigitsLength)
:m_Digits(DigitsSize(DigitsLength.GetDigits()),0)
{
m_Exponent=0;
m_Sign=0;
*this=DefultFloatValue;
}
TLargeFloat::TLargeFloat(long double FloatValue,const TCharacter&)//内部使用 浮点数转化为 TLargeFloat,并采用默认精度
:m_Digits(DigitsSize(emLongDoubleDigits),0)
{
m_Exponent=0;
m_Sign=0;
*this=FloatValue;
}
TLargeFloat::TLargeFloat(const SelfType& Value)
:m_Digits(Value.m_Digits)
{
m_Exponent=Value.m_Exponent;
m_Sign=Value.m_Sign;
}
TLargeFloat::SelfType& TLargeFloat::operator = (const SelfType& Value)
{
m_Digits=Value.m_Digits;
m_Exponent=Value.m_Exponent;
m_Sign=Value.m_Sign;
return *this;
}
void TLargeFloat::SetDigitsLength(unsigned int uiDigitsLength)//重新设置10进制有效位数
{
m_Digits.resize(DigitsSize(uiDigitsLength),0);
}
unsigned int TLargeFloat::GetDigitsLength() const//返回当前的10进制有效位数
{
return m_Digits.size()*em10Power;
}
void TLargeFloat::Clear()
{
//清零
int size=m_Digits.size();
for (int i=0;i<size;++i)
m_Digits[i]=0;
m_Exponent=0;
m_Sign=0;
}
TLargeFloat::SelfType& TLargeFloat::operator = (long double fValue)
{
Clear();
if (0==fValue)
{
//do nothing;
}
else
{
if (fValue>0)
m_Sign=1;
else
{
m_Sign=-1;
fValue=-fValue;
}
if (FloatIsInteger(fValue))//对 为"可表示整数" 的浮点数 进行特殊处理 无误差转换
{
long double tf=fValue;
int n=0;
for (;;++n)
{
if (0==tf) break;
tf/=emBase;
tf=floor(tf);
}
m_Exponent=n*em10Power;
for (int i=0;i<n;++i)
{
m_Digits[n-1-i]=Int32bit(TMaxInt(fValue)%emBase);
fValue=floor(fValue/emBase);
}
}
else//一般的浮点数 转化中可能产生小的误差
{
m_Exponent=int(floor(log10(fValue)))+1;//得到10为底的指数
fValue/=pow(10.0,(long double)(m_Exponent.AsInt()));
int size=m_Digits.size();
int minsize=Private_::min(size,emLongDoubleDigits/em10Power+1);
for (int i=0;i<minsize;++i)//得到小数位
{
if (0==fValue) break;
fValue*=emBase;
Int32bit IValue=Int32bit(floor(fValue));
fValue-=IValue;
m_Digits[i]=IValue;
if (i==minsize-1)
{
if (fValue*emBase*2>=emBase)//四舍五入
{
m_Digits[i]+=1;
for (int j=i;j>0;--j)
{
if (m_Digits[j]>=emBase)//进位
{
m_Digits[j-1]+=1;
m_Digits[j]-=emBase;
}
else
break;
}//for j
}//if
}//if
}
}//for i
}
Canonicity();
return *this;
}
long double TLargeFloat::AsFloat() const
{
//
if ( ((m_Exponent.AsInt())>=emLongDoubleMaxExponent)
||((m_Exponent.AsInt())<=emLongDoubleMinExponent) )
{
throw TException("ERROR:TLargeFloat::AsFloat()");
}
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);
}
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;
}
/////////////////////
