Hundredsoft

VC++で拡張倍精度浮動小数点形式をどうしても使いたい

Visual Stdioでは拡張倍精度浮動小数(80bit)が使えません。
long double と書いても、64bit精度になります。
通常、仮数部53bit(暗黙Bit含め)で不足することはないでしょうが、
どうしても仮数部64bitが必要という場合のサンプルクラスです。
様々なコンストラクタ、四則演算・比較の他、FPUでサポートしている三角関数や平方根を入れました。

/* (c) Hundredsoft Corporation 2010 All rights reserved.

	Ldouble.h -  VC 拡張倍精度クラス

*/

#ifndef LDOUBLE_INCLUDED
#define LDOUBLE_INCLUDED

#include <memory.h>

typedef struct {
	unsigned char d[12];
} LONG_DOUBLE_STR, *LPLONG_DOUBLE_STR;

void Ldouble_init()
{
	unsigned short fcw=0x033f;	// round to nearest, 64 bits, exceptions disabled
	__asm {
		fldcw     [fcw]
	}
}

class Ldouble
{
public:
	Ldouble() {
		memset(&m_data, 0, sizeof(m_data));
	}
	virtual ~Ldouble(){}

	Ldouble(const Ldouble& v) {
		if (this != &v){
			memcpy(&m_data, &v.m_data, sizeof(m_data));
		}
	}
	Ldouble(const __int64& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [v]
			fild qword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const unsigned __int64& v) {
		__int64 t = 0x8000000000000000;
		int i;
		for (i=63; i>=0; i--, t>>=1){
			if (v & t) break;
		}
		if (i < 0){
			memset(&m_data, 0, sizeof(m_data));
		}else{
			t = v;
			t <<= 63-i; // 仮数部
			int ex = i + 16383; // 指数部
			memcpy(&m_data, &t, sizeof(__int64));
			*(int*)&m_data.d[8] = ex;
		}
	}
	Ldouble(const int& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [v]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const unsigned int& v) {
		__int64 i = v;
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			lea  eax, [i]
			fild qword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const short& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		int i = v;
		__asm {
			lea  eax, [i]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const unsigned short& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		int i = v;
		__asm {
			lea  eax, [i]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const char& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		int i = v;
		__asm {
			lea  eax, [i]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const unsigned char& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		int i = v;
		__asm {
			lea  eax, [i]
			fild dword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble(const double& v) {
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [v]
			fld  qword ptr [eax]
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
	}
	Ldouble& operator=(const Ldouble& v){
		if (this != &v){
			memcpy(&m_data, &v.m_data, sizeof(m_data));
		}
		return *this;
	}

	// Unary -								
	Ldouble operator-(){
		Ldouble v = *this;
		v.m_data.d[9] ^= 0x80;
		return v;
	}
	Ldouble operator++(){
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fadd
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble operator++(int){
		Ldouble v = *this;
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fadd
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return v;
	}
	Ldouble operator--(){
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fsub
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble operator--(int){
		Ldouble v = *this;
		LPLONG_DOUBLE_STR m = &m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fsub
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return v;
	}

	// Binary operator +=, -=, *=, /=
	Ldouble& operator+=(const Ldouble& v){
		LPLONG_DOUBLE_STR m = &m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			fadd
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble& operator-=(const Ldouble& v){
		LPLONG_DOUBLE_STR m = &m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			fsub
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble& operator*=(const Ldouble& v){
		LPLONG_DOUBLE_STR m = &m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			fmul
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}
	Ldouble& operator/=(const Ldouble& v){
		LPLONG_DOUBLE_STR m = &m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			fdiv
			mov  eax, [m]
			fstp tbyte ptr [eax]
		}
		return *this;
	}

	// compare
	// this <  v: -
	// this == v: 0
	// this >  v: +
	int comp(const Ldouble& v) const{
		LPLONG_DOUBLE_STR m = (LPLONG_DOUBLE_STR)&m_data;
		LPLONG_DOUBLE_STR mv = (LPLONG_DOUBLE_STR)&v.m_data;
		short sts;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			mov  eax, [mv]
			fld  tbyte ptr [eax]
			lea  eax, [sts]
			fcompp
			fnstsw word ptr [eax]
		}
		return ((sts & 0x4100) - 1) ^ 0x3fff;
	}
	bool operator==(const Ldouble& v) const{
		return (comp(v) == 0);
	}
	bool operator!=(const Ldouble& v) const{
		return (comp(v) != 0);
	}
	bool operator<(const Ldouble& v) const{
		return (comp(v) < 0);
	}
	bool operator<=(const Ldouble& v) const{
		return (comp(v) <= 0);
	}
	bool operator>(const Ldouble& v) const{
		return (comp(v) > 0);
	}
	bool operator>=(const Ldouble& v) const{
		return (comp(v) >= 0);
	}

	// cast
	operator double() {
		LPLONG_DOUBLE_STR m = (LPLONG_DOUBLE_STR)&m_data;
		double t;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fstp qword ptr [t]
		}
		return t;
	}
	operator int() {
		LPLONG_DOUBLE_STR m = (LPLONG_DOUBLE_STR)&m_data;
		double t;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fstp qword ptr [t]
		}
		return (int)t;
	}

	//friend Binary operator	 +,-,*,/
	friend  const Ldouble operator+(const  Ldouble&a, const Ldouble&b){
		Ldouble p(a);
		p += b;
		return p;
	}
	friend  const Ldouble operator-(const  Ldouble&a, const Ldouble&b){
		Ldouble p(a);
		p -= b;
		return p;
	}
	friend  const Ldouble operator*(const  Ldouble&a, const Ldouble&b){
		Ldouble p(a);
		p *= b;
		return p;
	}
	friend  const Ldouble operator/(const  Ldouble&a, const Ldouble&b){
		Ldouble p(a);
		p /= b;
		return p;
	}


	LONG_DOUBLE_STR	m_data;
};

namespace LD {
	Ldouble abs(const Ldouble& s){
		Ldouble v(s);
		v.m_data.d[9] &= 0x7f;
		return v;
	}

	Ldouble sin(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fsin
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble cos(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fcos
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble tan(const Ldouble& s){
		Ldouble v(s);
		double t;
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fptan
			lea  ebx, [t]
			fstp qword ptr [ebx]
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble atan(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fld1
			fpatan
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble atan2(const Ldouble& x, const Ldouble& y){
		Ldouble v(y);
		LPLONG_DOUBLE_STR xm = (LPLONG_DOUBLE_STR)&x.m_data;
		LPLONG_DOUBLE_STR vm = &v.m_data;
		__asm {
			mov  eax, [xm]
			fld  tbyte ptr [eax]
			mov  eax, [vm]
			fld  tbyte ptr [eax]
			fpatan
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble sqrt(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fsqrt
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble pi(){
		Ldouble v;
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			mov  eax, [m]
			fldpi
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble log(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			fldln2
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fyl2x
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble log10(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		__asm {
			fldlg2
			mov  eax, [m]
			fld  tbyte ptr [eax]
			fyl2x
			fstp tbyte ptr [eax]
		}
		return v;
	}

	Ldouble exp(const Ldouble& s){
		Ldouble v(s);
		LPLONG_DOUBLE_STR m = &v.m_data;
		double half = 0.5;
		bool sign = ((m->d[9] & 0x80) != 0);
		m->d[9] &= 0x7f;

		__asm {
			fldl2e                 // log2e
			mov  eax, [m]
			fld  tbyte ptr [eax]   // log2e, v
			fmul                   // log2e*v
			fld  st(0)             // log2e*v, log2e*v
			lea  ebx, [half]
			fld  qword ptr [ebx]   // log2e*v, log2e*v, 0.5
			fsub                   // log2e*v, log2e*v-0.5
			frndint                // log2e*v, (int)(log2e*v-0.5)
			fxch                   // (int)(log2e*v-0.5), log2e*v
			fsub  st(0), st(1)     // (int)(log2e*v-0.5), log2e*v - (int)(log2e*v-0.5)
			f2xm1                  // (int)(log2e*v-0.5), 2^(log2e*v - (int)(log2e*v-0.5))-1
			fld1                   // (int)(log2e*v-0.5), 2^(log2e*v - (int)(log2e*v-0.5))-1, 1
			fadd                   // (int)(log2e*v-0.5), 2^(log2e*v - (int)(log2e*v-0.5))
			fscale                 // (int)(log2e*v-0.5), 2^(log2e*v - (int)(log2e*v-0.5)) * 2^(int)(log2e*v-0.5)
			fstp tbyte ptr [eax]   // (int)(log2e*v-0.5)
			fstp qword ptr [ebx]   // 
		}

		if (sign){
			__asm {
				fld1
				mov  eax, [m]
				fld  tbyte ptr [eax]
				fdiv
				fstp tbyte ptr [eax]
			}
		}
		return v;
	}

	Ldouble pow(const Ldouble& x, const Ldouble& y){
		Ldouble u(x),v(y);
		LPLONG_DOUBLE_STR mx = &u.m_data;
		LPLONG_DOUBLE_STR m = &v.m_data;
		double half = 0.5;
		bool xsign = ((mx->d[9] & 0x80) != 0);
		bool ysign = ((m->d[9] & 0x80) != 0);
		if (xsign){
			int ny = (int)v;
			Ldouble t = ny;
			if (t != y){
				throw "Invalid value.";
			}
			xsign = (ny & 1) ? true : false;
		}
		mx->d[9] &= 0x7f;
		m->d[9] &= 0x7f;

		__asm {
			fld1                   // 1
			mov  eax, [mx]
			fld  tbyte ptr [eax]   // 1, x
			fyl2x                  // log2x
			mov  eax, [m]
			fld  tbyte ptr [eax]   // log2x, v
			fmul                   // log2x*v
			fld  st(0)             // log2x*v, log2x*v
			lea  ebx, [half]
			fld  qword ptr [ebx]   // log2x*v, log2x*v, 0.5
			fsub                   // log2x*v, log2x*v-0.5
			frndint                // log2x*v, (int)(log2x*v-0.5)
			fxch                   // (int)(log2x*v-0.5), log2x*v
			fsub  st(0), st(1)     // (int)(log2x*v-0.5), log2x*v - (int)(log2x*v-0.5)
			f2xm1                  // (int)(log2x*v-0.5), 2^(log2x*v - (int)(log2x*v-0.5))-1
			fld1                   // (int)(log2x*v-0.5), 2^(log2x*v - (int)(log2x*v-0.5))-1, 1
			fadd                   // (int)(log2x*v-0.5), 2^(log2x*v - (int)(log2x*v-0.5))
			fscale                 // (int)(log2x*v-0.5), 2^(log2x*v - (int)(log2x*v-0.5)) * 2^(int)(log2x*v-0.5)
			fstp tbyte ptr [eax]   // (int)(log2x*v-0.5)
			fstp qword ptr [ebx]   // 
		}

		if (ysign){
			__asm {
				fld1
				mov  eax, [m]
				fld  tbyte ptr [eax]
				fdiv
				fstp tbyte ptr [eax]
			}
		}
		if (xsign){
			m->d[9] |= 0x80;
		}
		return v;
	}
} // end of namespace

#endif

Ldoubleクラスを使う前に、Ldouble_init()を呼び出して、FPUに拡張倍精度を指示する必要があります。
このコードは動作しますが、あまり速くありません。
実行速度についてはコンパイラが拡張倍精度をサポートしてくれないと、どうにもなりません。
実行速度をそれほど求めていないが、仮数部64bitの精度が必要であるか、
最大値が10進で300桁以上になる場合（拡張倍精度ではおよそ5000桁まで可能）にのみ有効です。
Win32アプリで拡張倍精度計算を行い、同時に実行速度も求めるなら、VC以外のコンパイラ(C++ Builder等)を利用した方が良いでしょう。

---

Tags: プログラムメモ