分享
 
 
 

自己动手打造“超高精度浮点数类”

王朝other·作者佚名  2007-01-28
窄屏简体版  字體: |||超大  

自己动手打造“超高精度浮点数类”

HouSisong@263.net

很多人可能都想自己写一个能够执行任意精度计算的浮点数;: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;

}

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

 
 
 
免责声明:本文为网络用户发布,其观点仅代表作者个人观点,与本站无关,本站仅提供信息存储服务。文中陈述内容未经本站证实,其真实性、完整性、及时性本站不作任何保证或承诺,请读者仅作参考,并请自行核实相关内容。
2023年上半年GDP全球前十五强
 百态   2023-10-24
美众议院议长启动对拜登的弹劾调查
 百态   2023-09-13
上海、济南、武汉等多地出现不明坠落物
 探索   2023-09-06
印度或要将国名改为“巴拉特”
 百态   2023-09-06
男子为女友送行,买票不登机被捕
 百态   2023-08-20
手机地震预警功能怎么开?
 干货   2023-08-06
女子4年卖2套房花700多万做美容:不但没变美脸,面部还出现变形
 百态   2023-08-04
住户一楼被水淹 还冲来8头猪
 百态   2023-07-31
女子体内爬出大量瓜子状活虫
 百态   2023-07-25
地球连续35年收到神秘规律性信号,网友:不要回答!
 探索   2023-07-21
全球镓价格本周大涨27%
 探索   2023-07-09
钱都流向了那些不缺钱的人,苦都留给了能吃苦的人
 探索   2023-07-02
倩女手游刀客魅者强控制(强混乱强眩晕强睡眠)和对应控制抗性的关系
 百态   2020-08-20
美国5月9日最新疫情:美国确诊人数突破131万
 百态   2020-05-09
荷兰政府宣布将集体辞职
 干货   2020-04-30
倩女幽魂手游师徒任务情义春秋猜成语答案逍遥观:鹏程万里
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案神机营:射石饮羽
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案昆仑山:拔刀相助
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案天工阁:鬼斧神工
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案丝路古道:单枪匹马
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:与虎谋皮
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:李代桃僵
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:指鹿为马
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案金陵:小鸟依人
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案金陵:千金买邻
 干货   2019-11-12
 
推荐阅读
 
 
 
>>返回首頁<<
 
靜靜地坐在廢墟上,四周的荒凉一望無際,忽然覺得,淒涼也很美
© 2005- 王朝網路 版權所有