📄 高精度浮点数.c
字号:
类的声明文件: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()");
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -