mp: add mptod and dtomp

3 files changed, 105 insertions, 3 deletions

diff --git a/sys/man/2/mp b/sys/man/2/mp
index 678f0eb7f..ed757d8f9 100644
--- a/sys/man/2/mp
+++ b/sys/man/2/mp
@@ -1,6 +1,6 @@
 .TH MP 2
 .SH NAME
-mpsetminbits, mpnew, mpfree, mpbits, mpnorm, mpcopy, mpassign, mprand, mpnrand, strtomp, mpfmt,mptoa, betomp, mptobe, mptober, letomp, mptole, mptolel, mptoui, uitomp, mptoi, itomp, uvtomp, mptouv, vtomp, mptov, mpdigdiv, mpadd, mpsub, mpleft, mpright, mpmul, mpexp, mpmod, mpmodadd, mpmodsub, mpmodmul, mpdiv, mpcmp, mpsel, mpextendedgcd, mpinvert, mpsignif, mplowbits0, mpvecdigmuladd, mpvecdigmulsub, mpvecadd, mpvecsub, mpveccmp, mpvecmul, mpmagcmp, mpmagadd, mpmagsub, crtpre, crtin, crtout, crtprefree, crtresfree \- extended precision arithmetic
+mpsetminbits, mpnew, mpfree, mpbits, mpnorm, mpcopy, mpassign, mprand, mpnrand, strtomp, mpfmt,mptoa, betomp, mptobe, mptober, letomp, mptole, mptolel, mptoui, uitomp, mptoi, itomp, uvtomp, mptouv, vtomp, mptov, mptod, dtomp, mpdigdiv, mpadd, mpsub, mpleft, mpright, mpmul, mpexp, mpmod, mpmodadd, mpmodsub, mpmodmul, mpdiv, mpcmp, mpsel, mpextendedgcd, mpinvert, mpsignif, mplowbits0, mpvecdigmuladd, mpvecdigmulsub, mpvecadd, mpvecsub, mpveccmp, mpvecmul, mpmagcmp, mpmagadd, mpmagsub, crtpre, crtin, crtout, crtprefree, crtresfree \- extended precision arithmetic
 .SH SYNOPSIS
 .B #include <u.h>
 .br
@@ -88,6 +88,12 @@ mpint*	uvtomp(uvlong, mpint*)
 uvlong	mptouv(mpint*)
 .PP
 .B
+mpint*	dtomp(double, mpint*)
+.PP
+.B
+double	mptod(mpint*)
+.PP
+.B
 void	mpadd(mpint *b1, mpint *b2, mpint *sum)
 .PP
 .B
@@ -271,8 +277,9 @@ This includes
 .IR strtomp ,
 .IR itomp ,
 .IR uitomp ,
+.IR btomp ,
 and
-.IR btomp .
+.IR dtomp .
 These functions, in addition to
 .I mpnew
 and
@@ -474,7 +481,7 @@ is
 .BR nil ,
 a new integer is allocated and returned as the result.
 .PP
-The integer conversions are:
+The integer (and floating point) conversions are:
 .TF Mptouv
 .TP
 .I mptoui
@@ -500,12 +507,23 @@ The integer conversions are:
 .TP
 .I vtomp
 .BR "vlong" -> mpint
+.TP
+.I mptod
+.BR mpint -> "double"
+.TP
+.I dtomp
+.BR "double" -> mpint
 .PD
 .PP
 When converting to the base integer types, if the integer is too large,
 the largest integer of the appropriate sign
 and size is returned.
 .PP
+When converting to and from floating point, results are rounded using IEEE 754 "round to nearest".
+If the integer is too large in magnitude,
+.I mptod 
+returns infinity of the appropriate sign.
+.PP
 The mathematical functions are:
 .TF mpmagadd
 .TP
diff --git a/sys/src/libmp/port/mkfile b/sys/src/libmp/port/mkfile
index 9aaad2d80..a94c4c92d 100644
--- a/sys/src/libmp/port/mkfile
+++ b/sys/src/libmp/port/mkfile
@@ -42,6 +42,7 @@ FILES=\
 	cnfield\
 	gmfield\
 	mplogic\
+	mptod\
 
 ALLOFILES=${FILES:%=%.$O}
 # cull things in the per-machine directories from this list
diff --git a/sys/src/libmp/port/mptod.c b/sys/src/libmp/port/mptod.c
new file mode 100644
index 000000000..d74e702fb
--- /dev/null
+++ b/sys/src/libmp/port/mptod.c
@@ -0,0 +1,83 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+double
+mptod(mpint *a)
+{
+	u64int v;
+	mpdigit w, r;
+	int sf, i, n, m, s;
+	FPdbleword x;
+	
+	if(a->top == 0) return 0.0;
+	sf = mpsignif(a);
+	if(sf > 1024) return Inf(a->sign);
+	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 & (1<<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 & (1<<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 Inf(a->sign);
+	}
+	x.lo = v;
+	x.hi = (u32int)(v >> 32) & (1<<20) - 1 | sf + 1022 << 20 | a->sign & 1<<31;
+	return x.x;
+}
+
+mpint *
+dtomp(double d, mpint *a)
+{
+	FPdbleword x;
+	uvlong 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 | (uvlong)(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;
+}