shithub: sl

ref: e7ec1e13f5e1814f647bcbc18e73eb6be5ce4f16
dir: /3rd/mp/mptod.c/

#include "platform.h"
#include "mp.h"

extern double D_PINF, D_NINF;

double
mptod(mpint *a)
{
	uint64_t v;
	mpdigit w, r;
	int sf, i, n;
	uint32_t 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;
			v = 0;
		}else{
			v <<= m;
		}
		s = (Dbits - m) & (Dbits - 1);
		v |= 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;
}