shithub: riscv

Info  •  Files  •  Log  •  Branches

mp: add mptod and dtomp

--- 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 @@

 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 @@

 .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 @@

 .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,11 +507,22 @@

 .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

--- a/sys/src/libmp/port/mkfile

+++ b/sys/src/libmp/port/mkfile

@@ -42,6 +42,7 @@

 	cnfield\

 	gmfield\

 	mplogic\

+	mptod\

 ALLOFILES=${FILES:%=%.$O}

 # cull things in the per-machine directories from this list

--- /dev/null

+++ b/sys/src/libmp/port/mptod.c

@@ -1,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;

+}