summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/get_ld.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/get_ld.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/get_ld.c227
1 files changed, 142 insertions, 85 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/get_ld.c b/Build/source/libs/mpfr/mpfr-src/src/get_ld.c
index b1d77c77145..8a3ad88ad57 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/get_ld.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/get_ld.c
@@ -21,94 +21,28 @@ along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
-#include <float.h>
+#include <float.h> /* needed so that MPFR_LDBL_MANT_DIG is correctly defined */
#include "mpfr-impl.h"
-#ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE
+#if defined(HAVE_LDOUBLE_IS_DOUBLE)
+/* special code when "long double" is the same format as "double" */
long double
mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
-
- if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
- return (long double) mpfr_get_d (x, rnd_mode);
- else /* now x is a normal non-zero number */
- {
- long double r; /* result */
- long double m;
- double s; /* part of result */
- mpfr_exp_t sh; /* exponent shift, so that x/2^sh is in the double range */
- mpfr_t y, z;
- int sign;
-
- /* first round x to the target long double precision, so that
- all subsequent operations are exact (this avoids double rounding
- problems) */
- mpfr_init2 (y, MPFR_LDBL_MANT_DIG);
- mpfr_init2 (z, MPFR_LDBL_MANT_DIG);
- /* Note about the precision of z: even though IEEE_DBL_MANT_DIG is
- sufficient, z has been set to the same precision as y so that
- the mpfr_sub below calls mpfr_sub1sp, which is faster than the
- generic subtraction, even in this particular case (from tests
- done by Patrick Pelissier on a 64-bit Core2 Duo against r7285).
- But here there is an important cancellation in the subtraction.
- TODO: get more information about what has been tested. */
-
- mpfr_set (y, x, rnd_mode);
- sh = MPFR_GET_EXP (y);
- sign = MPFR_SIGN (y);
- MPFR_SET_EXP (y, 0);
- MPFR_SET_POS (y);
-
- r = 0.0;
- do {
- s = mpfr_get_d (y, MPFR_RNDN); /* high part of y */
- r += (long double) s;
- mpfr_set_d (z, s, MPFR_RNDN); /* exact */
- mpfr_sub (y, y, z, MPFR_RNDN); /* exact */
- } while (!MPFR_IS_ZERO (y));
-
- mpfr_clear (z);
- mpfr_clear (y);
-
- /* we now have to multiply back by 2^sh */
- MPFR_ASSERTD (r > 0);
- if (sh != 0)
- {
- /* An overflow may occurs (example: 0.5*2^1024) */
- while (r < 1.0)
- {
- r += r;
- sh--;
- }
-
- if (sh > 0)
- m = 2.0;
- else
- {
- m = 0.5;
- sh = -sh;
- }
-
- for (;;)
- {
- if (sh % 2)
- r = r * m;
- sh >>= 1;
- if (sh == 0)
- break;
- m = m * m;
- }
- }
- if (sign < 0)
- r = -r;
- return r;
- }
+ return (long double) mpfr_get_d (x, rnd_mode);
}
-#else
+#elif defined(HAVE_LDOUBLE_IEEE_EXT_LITTLE)
+/* Note: The code will return a result with a 64-bit precision, even
+ if the rounding precision is only 53 bits like on FreeBSD and
+ NetBSD 6- (or with GCC's -mpc64 option to simulate this on other
+ platforms). This is consistent with how strtold behaves in these
+ cases, for instance. */
+
+/* special code for IEEE 754 little-endian extended format */
long double
mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
{
@@ -122,7 +56,7 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG);
inex = mpfr_set (tmp, x, rnd_mode);
- mpfr_set_emin (-16381-63); /* emin=-16444 */
+ mpfr_set_emin (-16381-63); /* emin=-16444, see below */
mpfr_set_emax (16384);
mpfr_subnormalize (tmp, mpfr_check_range (tmp, inex, rnd_mode), rnd_mode);
mpfr_prec_round (tmp, 64, MPFR_RNDZ); /* exact */
@@ -135,10 +69,12 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
tmpmant = MPFR_MANT (tmp);
e = MPFR_GET_EXP (tmp);
- /* the smallest normal number is 2^(-16382), which is 0.5*2^(-16381)
- in MPFR, thus any exponent <= -16382 corresponds to a subnormal
- number */
+ /* The smallest positive normal number is 2^(-16382), which is
+ 0.5*2^(-16381) in MPFR, thus any exponent <= -16382 corresponds to a
+ subnormal number. The smallest positive subnormal number is 2^(-16445)
+ which is 0.5*2^(-16444) in MPFR thus 0 <= denorm <= 63. */
denorm = MPFR_UNLIKELY (e <= -16382) ? - e - 16382 + 1 : 0;
+ MPFR_ASSERTD (0 <= denorm && denorm < 64);
#if GMP_NUMB_BITS >= 64
ld.s.manl = (tmpmant[0] >> denorm);
ld.s.manh = (tmpmant[0] >> denorm) >> 32;
@@ -153,7 +89,7 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
ld.s.manl = (tmpmant[0] >> denorm) | (tmpmant[1] << (32 - denorm));
ld.s.manh = tmpmant[1] >> denorm;
}
- else /* 32 <= denorm <= 64 */
+ else /* 32 <= denorm < 64 */
{
ld.s.manl = tmpmant[1] >> (denorm - 32);
ld.s.manh = 0;
@@ -177,6 +113,128 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
return ld.ld;
}
+#else
+
+/* generic code */
+long double
+mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
+{
+ if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+ return (long double) mpfr_get_d (x, rnd_mode);
+ else /* now x is a normal non-zero number */
+ {
+ long double r; /* result */
+ double s; /* part of result */
+ MPFR_SAVE_EXPO_DECL (expo);
+
+ MPFR_SAVE_EXPO_MARK (expo);
+
+#if defined(HAVE_LDOUBLE_MAYBE_DOUBLE_DOUBLE)
+ if (MPFR_LDBL_MANT_DIG == 106)
+ {
+ /* Assume double-double format (as found with the PowerPC ABI).
+ The generic code below isn't used because numbers with
+ precision > 106 would not be supported. */
+ s = mpfr_get_d (x, MPFR_RNDN); /* high part of x */
+ /* Let's first consider special cases separately. The test for
+ infinity is really needed to avoid a NaN result. The test
+ for NaN is mainly for optimization. The test for 0 is useful
+ to get the correct sign (assuming mpfr_get_d supports signed
+ zeros on the implementation). */
+ if (s == 0 || DOUBLE_ISNAN (s) || DOUBLE_ISINF (s))
+ r = (long double) s;
+ else
+ {
+ mpfr_t y, z;
+
+ mpfr_init2 (y, mpfr_get_prec (x));
+ mpfr_init2 (z, IEEE_DBL_MANT_DIG); /* keep the precision small */
+ mpfr_set_d (z, s, MPFR_RNDN); /* exact */
+ mpfr_sub (y, x, z, MPFR_RNDN); /* exact */
+ /* Add the second part of y (in the correct rounding mode). */
+ r = (long double) s + (long double) mpfr_get_d (y, rnd_mode);
+ mpfr_clear (z);
+ mpfr_clear (y);
+ }
+ }
+ else
+#endif
+ {
+ long double m;
+ mpfr_exp_t sh; /* exponent shift -> x/2^sh is in the double range */
+ mpfr_t y, z;
+ int sign;
+
+ /* First round x to the target long double precision, so that
+ all subsequent operations are exact (this avoids double rounding
+ problems). However if the format contains numbers that have more
+ precision, MPFR won't be able to generate such numbers. */
+ mpfr_init2 (y, MPFR_LDBL_MANT_DIG);
+ mpfr_init2 (z, MPFR_LDBL_MANT_DIG);
+ /* Note about the precision of z: even though IEEE_DBL_MANT_DIG is
+ sufficient, z has been set to the same precision as y so that
+ the mpfr_sub below calls mpfr_sub1sp, which is faster than the
+ generic subtraction, even in this particular case (from tests
+ done by Patrick Pelissier on a 64-bit Core2 Duo against r7285).
+ But here there is an important cancellation in the subtraction.
+ TODO: get more information about what has been tested. */
+
+ mpfr_set (y, x, rnd_mode);
+ sh = MPFR_GET_EXP (y);
+ sign = MPFR_SIGN (y);
+ MPFR_SET_EXP (y, 0);
+ MPFR_SET_POS (y);
+
+ r = 0.0;
+ do
+ {
+ s = mpfr_get_d (y, MPFR_RNDN); /* high part of y */
+ r += (long double) s;
+ mpfr_set_d (z, s, MPFR_RNDN); /* exact */
+ mpfr_sub (y, y, z, MPFR_RNDN); /* exact */
+ }
+ while (!MPFR_IS_ZERO (y));
+
+ mpfr_clear (z);
+ mpfr_clear (y);
+
+ /* we now have to multiply back by 2^sh */
+ MPFR_ASSERTD (r > 0);
+ if (sh != 0)
+ {
+ /* An overflow may occur (example: 0.5*2^1024) */
+ while (r < 1.0)
+ {
+ r += r;
+ sh--;
+ }
+
+ if (sh > 0)
+ m = 2.0;
+ else
+ {
+ m = 0.5;
+ sh = -sh;
+ }
+
+ for (;;)
+ {
+ if (sh % 2)
+ r = r * m;
+ sh >>= 1;
+ if (sh == 0)
+ break;
+ m = m * m;
+ }
+ }
+ if (sign < 0)
+ r = -r;
+ }
+ MPFR_SAVE_EXPO_FREE (expo);
+ return r;
+ }
+}
+
#endif
/* contributed by Damien Stehle */
@@ -190,8 +248,7 @@ mpfr_get_ld_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
return (long double) mpfr_get_d_2exp (expptr, src, rnd_mode);
- tmp[0] = *src; /* Hack copy mpfr_t */
- MPFR_SET_EXP (tmp, 0);
+ MPFR_ALIAS (tmp, src, MPFR_SIGN (src), 0);
ret = mpfr_get_ld (tmp, rnd_mode);
if (MPFR_IS_PURE_FP(src))