shithub: femtolisp

ref: 20488fc4915662e0ec50bb8c4ebc2609d6e110b4
dir: /3rd/mp/mptod.c/

#include "platform.h"

extern double D_PINF, D_NINF;

double
mptod(mpint *a)
{
	uint64_t v;
	mpdigit w, r;
	int sf, i, n, m, s;
	FPdbleword x;

	if(a->top == 0) return 0.0;
	sf = (int)mpsignif(a);
	if(sf > 1024) return a->sign < 0 ? D_NINF : D_PINF;
	i = a->top - 1;
	v = a->p[i];
	n = sf & Dbits - 1;
	n |= n - 1 & Dbits;
	r = 0;
	if(n > 54){
		s = n - 54;
		r = v & (1ULL<<s) - 1;
		v >>= s;
	}
	while(n < 54){
		if(--i < 0)
			w = 0;
		else
			w = a->p[i];
		m = 54 - n;
		if(m > Dbits) m = Dbits;
		s = Dbits - m & Dbits - 1;
		v = v << m | w >> s;
		r = w & (1ULL<<s) - 1;
		n += m;
	}
	if((v & 3) == 1){
		while(--i >= 0)
			r |= a->p[i];
		if(r != 0)
			v++;
	}else
		v++;
	v >>= 1;
	while((v >> 53) != 0){
		v >>= 1;
		if(++sf > 1024)
			return a->sign < 0 ? D_NINF : D_PINF;
	}
	x.lo = v;
	x.hi = (uint32_t)(v >> 32) & (1<<20) - 1 | (sf + 1022) << 20 | a->sign & 1U<<31;
	return x.x;
}

mpint *
dtomp(double d, mpint *a)
{
	FPdbleword x;
	uint64_t v;
	int e;

	if(a == nil)
		a = mpnew(0);
	x.x = d;
	e = x.hi >> 20 & 2047;
	assert(e != 2047);
	if(e < 1022){
		mpassign(mpzero, a);
		return a;
	}
	v = x.lo | (uint64_t)(x.hi & (1<<20) - 1) << 32 | 1ULL<<52;
	if(e < 1075){
		v += (1ULL<<(1074 - e)) - (~v >> (1075 - e) & 1);
		v >>= 1075 - e;
	}
	uvtomp(v, a);
	if(e > 1075)
		mpleft(a, e - 1075, a);
	if((int)x.hi < 0)
		a->sign = -1;
	return a;
}