diff options
| -rw-r--r-- | configure.in | 22 | ||||
| -rw-r--r-- | get_d64.c | 120 | ||||
| -rw-r--r-- | mpfr-impl.h | 3 | ||||
| -rw-r--r-- | mpfr.h | 16 | ||||
| -rw-r--r-- | set_d64.c | 48 | ||||
| -rw-r--r-- | tests/tget_set_d64.c | 18 |
6 files changed, 180 insertions, 47 deletions
diff --git a/configure.in b/configure.in index 06f5c0885..88337c35e 100644 --- a/configure.in +++ b/configure.in @@ -48,6 +48,7 @@ AC_ARG_WITH(gmp_build, GMP_CFLAGS=`grep -w "CFLAGS =" $withval/Makefile | sed 's/CFLAGS = //'` GMP_CC=`grep -w "CC =" $withval/Makefile | sed 's/CC = //'` fi + use_gmp_build=yes ]) AC_ARG_WITH(gmp, [ --with-gmp=DIR GMP install directory ], [ @@ -88,7 +89,26 @@ AC_ARG_ENABLE(thread-safe, AC_ARG_ENABLE(decimal-float, [ --enable-decimal-float build conversion functions from/to decimal floats [[default=no]]], [ case $enableval in - yes) AC_DEFINE([MPFR_WANT_DECIMAL_FLOATS],1,[Build decimal float functions]) ;; + yes) AC_DEFINE([MPFR_WANT_DECIMAL_FLOATS],1, + [Build decimal float functions]) + AC_MSG_CHECKING(if compiler knows _Decimal64) + AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[_Decimal64 x;]])], + [AC_MSG_RESULT(yes) + if test "$use_gmp_build" != yes ; then + AC_MSG_ERROR([decimal float support requires --with-gmp-build]) + fi + AC_MSG_CHECKING(if _GMP_IEEE_FLOATS is defined) + AC_COMPILE_IFELSE([AC_LANG_SOURCE([[ +#include "gmp.h" +#include "gmp-impl.h" +#ifndef _GMP_IEEE_FLOATS +#error "_GMP_IEEE_FLOATS is not defined" +#endif]])],[AC_MSG_RESULT(yes)],[AC_MSG_RESULT(no) + AC_MSG_ERROR([decimal float support requires _GMP_IEEE_FLOATS])]) + ], + [AC_MSG_ERROR([Compiler doesn't know _Decimal64; try GCC >= 4.2] + )]) + ;; no) ;; *) AC_MSG_ERROR([bad value for --enable-decimal-float: yes or no]) ;; esac]) @@ -27,43 +27,126 @@ MA 02110-1301, USA. */ #include <string.h> /* for strlen */ #include "mpfr-impl.h" -#if WANT_DECIMAL64 +#if MPFR_WANT_DECIMAL_FLOATS + +#ifdef DPD_FORMAT +static _Decimal64 +bid_to_dpd (_Decimal64 d) +{ + union ieee_double_extract x; + union ieee_double_decimal64 y; + mp_limb_t h, l; + unsigned long d0, d1, d2, d3, d4, d5; /* declets: 0 <= d0..d5 < 999 */ + + y.d64 = d; + x.d = y.d; + G = x.s.exp << 2; /* x.s.exp is the biased exponent, in [0, 2048) */ +#if BITS_PER_MP_LIMB == 32 + G |= x.s.manh >> 18; /* manh has 20 bits */ +#elif BITS_PER_MP_LIMB == 64 + G |= x.s.manl >> 50; /* manl has 52 bits */ +#else +#error "wrong value of BITS_PER_MP_LIMB" +#endif + /* now G has 13 bits: 0 <= G < 8192 */ + Gh = G >> 8; + if (Gh >= 30) /* NaN or Inf have the same encoding in DPD or BID */ + return d; + + if (Gh >= 24) + { + h = 8 | (G & 1); /* 8 + G[w+4] */ + exp = (G >> 1) & 1023; + } + else /* 0 <= Gh < 24 */ + { + h = G & 7; + exp = G >> 3; + } + /* now significand is h*2^50 plus the remaining 50 bits */ +#if BITS_PER_MP_LIMB == 32 + h = h << 18 | (x.s.manh & 262143); + l = x.s.manl; + /* 2^32*h+l < 10^16 i.e. h <= 2328306 */ + d5 = (296 * h + l) % 1000; /* 2^32 = 296 mod 1000 */ + sub_ddmmss (h, l, h, l, 0, d5); + /* now h*2^32 + l is an exact multiple of 1000 */ + /* first divide by 8 */ + l = ((h & 7) << 29) | (l >> 3); + h = h >> 3; + /* now divide exactly by 125 */ + l = l * 652835029; /* 1/125 mod 2^32 */ + h = h / 125; + /* now 2^32*h+l < 10^13 i.e. h <= 2328 */ + d4 = (296 * h + l) % 1000; + sub_ddmmss (h, l, h, l, 0, d4); + l = ((h & 7) << 29) | (l >> 3); + h = h >> 3; + l = l * 652835029; + h = h / 125; + /* now 2^32*h+l < 10^10 i.e. h <= 2 */ + d3 = (296 * h + l) % 1000; + sub_ddmmss (h, l, h, l, 0, d4); + l = (h << 29) | (l >> 3); + l = l * 652835029; + /* now l < 10^7 */ +#else + l = h << 50 | (x.s.manl & 1125899906842623); + /* l < 10*2^50 theoretically, but l < 10^16 in practice */ + d5 = l % 1000; + l = l / 1000; /* l < 10^13 */ + d4 = l % 1000; + l = l / 1000; /* l < 10^10 */ + d3 = l % 1000; + l = l / 1000; /* l < 10^7 */ +#endif + d2 = l % 1000; + l = l / 1000; + d1 = l % 1000; + d0 = l / 10; + + /* now encode the declets */ +} +#endif /* DPD_FORMAT */ /* construct a decimal64 NaN */ -static decimal64 +static _Decimal64 get_decimal64_nan (void) { union ieee_double_extract x; + union ieee_double_decimal64 y; x.s.exp = 1984; /* G[0]..G[4] = 11111: quiet NaN */ - return x.d; + y.d = x.d; + return y.d64; } /* construct the decimal64 Inf with given sign */ -static decimal64 +static _Decimal64 get_decimal64_inf (int negative) { union ieee_double_extract x; + union ieee_double_decimal64 y; x.s.sig = (negative) ? 1 : 0; x.s.exp = 1920; /* G[0]..G[4] = 11110: Inf */ - return x.d; + y.d = x.d; + return y.d64; } /* construct the decimal64 zero with given sign */ -static decimal64 +static _Decimal64 get_decimal64_zero (int negative) { - union ieee_double_extract x; + union ieee_double_decimal64 y; - /* we assume here that +0 and -0 have the same encoding in - binary64 and decimal64 */ - x.d = negative ? DBL_NEG_ZERO : 0.0; - return x.d; + /* zero has the same representation in binary64 and decimal64 */ + y.d = negative ? DBL_NEG_ZERO : 0.0; + return y.d64; } /* construct the decimal64 smallest non-zero with given sign */ -static decimal64 +static _Decimal64 get_decimal64_min (int negative) { union ieee_double_extract x; @@ -78,7 +161,7 @@ get_decimal64_min (int negative) } /* construct the decimal64 largest finite number with given sign */ -static decimal64 +static _Decimal64 get_decimal64_max (int negative) { union ieee_double_extract x; @@ -94,10 +177,11 @@ get_decimal64_max (int negative) /* for BID format, the exponent exp is meant with a significand of the form [0.]sss...sss */ -static decimal64 +static _Decimal64 fill_decimal64 (int negative, mp_exp_t exp, char *s) { union ieee_double_extract x; + union ieee_double_decimal64 y; unsigned int i, l; mp_limb_t rp[2]; mp_size_t rn; @@ -135,15 +219,17 @@ fill_decimal64 (int negative, mp_exp_t exp, char *s) x.s.manl = (rp[0] ^ 9007199254740992) | ((exp & 1) << 51); #endif } - return x.d; + y.d = x.d; + return y.d64; } -decimal64 +_Decimal64 mpfr_get_decimal64 (mpfr_srcptr src, mp_rnd_t rnd_mode) { int negative; mp_exp_t e; + /* the encoding of NaN, Inf, zero is the same under DPD or BID */ if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src))) { if (MPFR_IS_NAN (src)) @@ -228,4 +314,4 @@ mpfr_get_decimal64 (mpfr_srcptr src, mp_rnd_t rnd_mode) } } -#endif /* WANT_DECIMAL64 */ +#endif /* MPFR_WANT_DECIMAL_FLOATS */ diff --git a/mpfr-impl.h b/mpfr-impl.h index fbebad7c5..0abe19675 100644 --- a/mpfr-impl.h +++ b/mpfr-impl.h @@ -355,7 +355,8 @@ typedef union ieee_double_extract Ieee_double_extract; # endif #endif - +/* to cast between binary64 and decimal64 */ +union ieee_double_decimal64 { double d; _Decimal64 d64; }; /****************************************************** *************** Long double macros ******************* @@ -124,14 +124,6 @@ typedef enum { MPFR_INF_KIND = 1, MPFR_ZERO_KIND = 2, MPFR_REGULAR_KIND = 3 } mpfr_kind_t; -/* the decimal64 conversion routines depend on _GMP_IEEE_FLOATS */ -#define WANT_DECIMAL64 (MPFR_WANT_DECIMAL_FLOATS && _GMP_IEEE_FLOATS) - -#if WANT_DECIMAL64 -/* FIXME: remove/replaced when decimal64 is a standard C type */ -typedef double decimal64; -#endif - /* GMP defines: + size_t: Standard size_t + __GMP_ATTRIBUTE_PURE Attribute for math functions. @@ -235,8 +227,8 @@ __MPFR_DECLSPEC void mpfr_set_default_prec _MPFR_PROTO((mpfr_prec_t)); __MPFR_DECLSPEC mp_prec_t mpfr_get_default_prec _MPFR_PROTO((void)); __MPFR_DECLSPEC int mpfr_set_d _MPFR_PROTO ((mpfr_ptr, double, mpfr_rnd_t)); -#if WANT_DECIMAL64 -__MPFR_DECLSPEC int mpfr_set_decimal64 _MPFR_PROTO ((mpfr_ptr, decimal64, +#if MPFR_WANT_DECIMAL_FLOATS +__MPFR_DECLSPEC int mpfr_set_decimal64 _MPFR_PROTO ((mpfr_ptr, _Decimal64, mp_rnd_t)); #endif __MPFR_DECLSPEC int @@ -292,8 +284,8 @@ __MPFR_DECLSPEC uintmax_t mpfr_get_uj _MPFR_PROTO ((mpfr_srcptr, mpfr_rnd_t)); __MPFR_DECLSPEC mp_exp_t mpfr_get_z_exp _MPFR_PROTO ((mpz_ptr, mpfr_srcptr)); __MPFR_DECLSPEC double mpfr_get_d _MPFR_PROTO ((mpfr_srcptr, mpfr_rnd_t)); -#if WANT_DECIMAL64 -__MPFR_DECLSPEC decimal64 mpfr_get_decimal64 _MPFR_PROTO ((mpfr_srcptr, +#if MPFR_WANT_DECIMAL_FLOATS +__MPFR_DECLSPEC _Decimal64 mpfr_get_decimal64 _MPFR_PROTO ((mpfr_srcptr, mp_rnd_t)); #endif __MPFR_DECLSPEC long double mpfr_get_ld _MPFR_PROTO ((mpfr_srcptr, @@ -28,11 +28,11 @@ MA 02110-1301, USA. */ #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" -#if WANT_DECIMAL64 +#if MPFR_WANT_DECIMAL_FLOATS #ifdef DPD_FORMAT -static decimal64 -dpd_to_bid (decimal64 d) +static _Decimal64 +dpd_to_bid (_Decimal64 d) { union ieee_double_extract x; unsigned int exp; @@ -125,7 +125,7 @@ dpd_to_bid (decimal64 d) if (Gh >= 24) { - d0 = 8 | (Gh & 2); /* 8 or 9 */ + d0 = 8 | (Gh & 1); /* 8 or 9 */ exp = (Gh >> 1) & 3; /* 0, 1, or 2 */ exp = (exp << 8) | (G & 255); /* 0 <= exp < 768 */ /* now convert exp, d0 to BID: the biased exponent is formed from @@ -175,10 +175,41 @@ dpd_to_bid (decimal64 d) } #endif /* DPD_FORMAT */ +#ifdef DEBUG +static void +print_decimal64 (_Decimal64 d) +{ + union ieee_double_extract x; + union ieee_double_decimal64 y; + unsigned int Gh, i; + + y.d64 = d; + x.d = y.d; + Gh = x.s.exp >> 6; + printf ("|%d|%d%d%d%d%d|", x.s.sig, Gh >> 4, (Gh >> 3) & 1, + (Gh >> 2) & 1, (Gh >> 1) & 1, Gh & 1); + printf ("%d%d%d%d%d%d|", (x.s.exp >> 5) & 1, (x.s.exp >> 4) & 1, + (x.s.exp >> 3) & 1, (x.s.exp >> 2) & 1, (x.s.exp >> 1) & 1, + x.s.exp & 1); +#if BITS_PER_MP_LIMB == 32 + for (i = 20; i > 0; i--) + printf ("%d", (x.s.manh >> (i - 1)) & 1); + printf ("|"); + for (i = 32; i > 0; i--) + printf ("%d", (x.s.manl >> (i - 1)) & 1); +#else + for (i = 52; i > 0; i--) + printf ("%d", (x.s.manl >> (i - 1)) & 1); +#endif + printf ("|\n"); +} +#endif + int -mpfr_set_decimal64 (mpfr_ptr r, decimal64 d, mp_rnd_t rnd_mode) +mpfr_set_decimal64 (mpfr_ptr r, _Decimal64 d, mp_rnd_t rnd_mode) { union ieee_double_extract x; + union ieee_double_decimal64 y; char s[23], *t; unsigned int Gh; /* most 5 significant bits from combination field */ int exp; /* exponent */ @@ -192,8 +223,9 @@ mpfr_set_decimal64 (mpfr_ptr r, decimal64 d, mp_rnd_t rnd_mode) 16 characters for mantissa, 1 character for E, 4 characters for exponent (including sign) */ - - x.d = d; + + y.d64 = d; + x.d = y.d; Gh = x.s.exp >> 6; if (Gh == 31) { @@ -266,4 +298,4 @@ mpfr_set_decimal64 (mpfr_ptr r, decimal64 d, mp_rnd_t rnd_mode) return mpfr_set_str (r, s, 10, rnd_mode); } -#endif /* WANT_DECIMAL64 */ +#endif /* MPFR_WANT_DECIMAL_FLOATS */ diff --git a/tests/tget_set_d64.c b/tests/tget_set_d64.c index f3fe98e59..6fd5d7331 100644 --- a/tests/tget_set_d64.c +++ b/tests/tget_set_d64.c @@ -22,14 +22,16 @@ MA 02110-1301, USA. */ #include <stdlib.h> /* for exit */ #include "mpfr-test.h" -#if WANT_DECIMAL64 +#if MPFR_WANT_DECIMAL_FLOATS static void -print_decimal64 (decimal64 d) +print_decimal64 (_Decimal64 d) { union ieee_double_extract x; + union ieee_double_decimal64 y; unsigned int Gh, i; - x.d = d; + y.d64 = d; + x.d = y.d; Gh = x.s.exp >> 6; printf ("|%d|%d%d%d%d%d|", x.s.sig, Gh >> 4, (Gh >> 3) & 1, (Gh >> 2) & 1, (Gh >> 1) & 1, Gh & 1); @@ -53,8 +55,8 @@ static void check_inf_nan () { mpfr_t x, y; - decimal64 d; - + _Decimal64 d; + mpfr_init2 (x, 123); mpfr_init2 (y, 123); @@ -170,7 +172,7 @@ static void check_random (void) { mpfr_t x, y; - decimal64 d; + _Decimal64 d; int i; mpfr_init2 (x, 49); @@ -197,7 +199,7 @@ check_random (void) mpfr_clear (x); mpfr_clear (y); } -#endif /* WANT_DECIMAL64 */ +#endif /* MPFR_WANT_DECIMAL_FLOATS */ int main (void) @@ -205,7 +207,7 @@ main (void) tests_start_mpfr (); mpfr_test_init (); -#if WANT_DECIMAL64 +#if MPFR_WANT_DECIMAL_FLOATS check_inf_nan (); check_random (); #endif |