<![CDATA[
/* (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
]]>
ソースはこちら