version 1.1.1.1, 2000/09/09 14:12:19 |
version 1.1.1.2, 2003/08/25 16:06:07 |
|
|
/* mpfr_set_d, mpfr_get_d -- convert a multiple precision floating-point number |
/* mpfr_set_d -- convert a machine double precision float to |
from/to a machine double precision float |
a multiple precision floating-point number |
|
|
Copyright (C) 1999 PolKA project, Inria Lorraine and Loria |
Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. |
|
|
This file is part of the MPFR Library. |
This file is part of the MPFR Library. |
|
|
The MPFR Library is free software; you can redistribute it and/or modify |
The MPFR Library is free software; you can redistribute it and/or modify |
it under the terms of the GNU Library General Public License as published by |
it under the terms of the GNU Lesser General Public License as published by |
the Free Software Foundation; either version 2 of the License, or (at your |
the Free Software Foundation; either version 2.1 of the License, or (at your |
option) any later version. |
option) any later version. |
|
|
The MPFR Library is distributed in the hope that it will be useful, but |
The MPFR Library is distributed in the hope that it will be useful, but |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |
License for more details. |
License for more details. |
|
|
You should have received a copy of the GNU Library General Public License |
You should have received a copy of the GNU Lesser General Public License |
along with the MPFR Library; see the file COPYING.LIB. If not, write to |
along with the MPFR Library; see the file COPYING.LIB. If not, write to |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
MA 02111-1307, USA. */ |
MA 02111-1307, USA. */ |
|
|
#if __GNUC__ /* gcc "patched" headers seem to omit isnan... */ |
|
extern int isnan(double); |
|
#endif |
|
#include <math.h> /* for isnan and NaN */ |
|
|
|
#include "gmp.h" |
#include "gmp.h" |
#include "gmp-impl.h" |
#include "gmp-impl.h" |
#include "longlong.h" |
#include "longlong.h" |
#include "mpfr.h" |
#include "mpfr.h" |
|
#include "mpfr-impl.h" |
|
|
#define NaN sqrt(-1) /* ensures a machine-independent NaN */ |
#if (BITS_PER_MP_LIMB==32) |
|
#define MPFR_LIMBS_PER_DOUBLE 2 |
|
#elif (BITS_PER_MP_LIMB >= 64) |
|
#define MPFR_LIMBS_PER_DOUBLE 1 |
|
#elif (BITS_PER_MP_LIMB == 16) |
|
#define MPFR_LIMBS_PER_DOUBLE 4 |
|
#endif |
|
|
|
static int __mpfr_extract_double _PROTO ((mp_ptr, double)); |
|
|
/* Included from gmp-2.0.2, patched to support denorms */ |
/* Included from gmp-2.0.2, patched to support denorms */ |
|
|
#ifdef XDEBUG |
#ifdef XDEBUG |
Line 42 extern int isnan(double); |
|
Line 46 extern int isnan(double); |
|
#define _GMP_IEEE_FLOATS 0 |
#define _GMP_IEEE_FLOATS 0 |
#endif |
#endif |
|
|
int |
static int |
#if __STDC__ |
__mpfr_extract_double (mp_ptr rp, double d) |
__mpfr_extract_double (mp_ptr rp, double d, int e) |
|
#else |
|
__mpfr_extract_double (rp, d, e) |
|
mp_ptr rp; |
|
double d; |
|
int e; |
|
#endif |
|
/* e=0 iff BITS_PER_MP_LIMB=32 and rp has only one limb */ |
/* e=0 iff BITS_PER_MP_LIMB=32 and rp has only one limb */ |
{ |
{ |
long exp; |
long exp; |
Line 70 __mpfr_extract_double (rp, d, e) |
|
Line 67 __mpfr_extract_double (rp, d, e) |
|
if (d == 0.0) |
if (d == 0.0) |
{ |
{ |
rp[0] = 0; |
rp[0] = 0; |
#if BITS_PER_MP_LIMB == 32 |
|
if (e) rp[1] = 0; |
|
#endif |
|
return 0; |
return 0; |
} |
} |
|
|
Line 82 __mpfr_extract_double (rp, d, e) |
|
Line 76 __mpfr_extract_double (rp, d, e) |
|
x.d = d; |
x.d = d; |
|
|
exp = x.s.exp; |
exp = x.s.exp; |
if (exp) |
if (exp) |
{ |
{ |
#if BITS_PER_MP_LIMB == 64 |
#if BITS_PER_MP_LIMB == 64 |
manl = (((mp_limb_t) 1 << 63) |
manl = ((MP_LIMB_T_ONE << 63) |
| ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11)); |
| ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11)); |
#else |
#else |
manh = ((mp_limb_t) 1 << 31) | (x.s.manh << 11) | (x.s.manl >> 21); |
manh = (MP_LIMB_T_ONE << 31) | (x.s.manh << 11) | (x.s.manl >> 21); |
manl = x.s.manl << 11; |
manl = x.s.manl << 11; |
#endif |
#endif |
} |
} |
else |
else /* denormalized number */ |
{ |
{ |
#if BITS_PER_MP_LIMB == 64 |
#if BITS_PER_MP_LIMB == 64 |
manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11); |
manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11); |
#else |
#else |
manh = (x.s.manh << 11) | (x.s.manl >> 21); |
manh = (x.s.manh << 11) /* high 21 bits */ |
manl = x.s.manl << 11; |
| (x.s.manl >> 21); /* middle 11 bits */ |
|
manl = x.s.manl << 11; /* low 21 bits */ |
#endif |
#endif |
} |
} |
} |
} |
Line 148 __mpfr_extract_double (rp, d, e) |
|
Line 143 __mpfr_extract_double (rp, d, e) |
|
} |
} |
#endif |
#endif |
|
|
if (exp) exp = (unsigned) exp - 1022; else exp = -1021; |
if (exp) exp = (unsigned) exp - 1022; else exp = -1021; |
|
|
#if BITS_PER_MP_LIMB == 64 |
#if BITS_PER_MP_LIMB == 64 |
rp[0] = manl; |
rp[0] = manl; |
#else |
#else |
if (e) { |
rp[1] = manh; |
rp[1] = manh; |
rp[0] = manl; |
rp[0] = manl; |
|
} |
|
else { |
|
rp[0] = manh; |
|
} |
|
#endif |
#endif |
|
|
return exp; |
return exp; |
} |
} |
|
|
/* End of part included from gmp-2.0.2 */ |
/* End of part included from gmp-2.0.2 */ |
/* Part included from gmp temporary releases */ |
|
double |
int |
#if __STDC__ |
mpfr_set_d (mpfr_ptr r, double d, mp_rnd_t rnd_mode) |
__mpfr_scale2 (double d, int exp) |
|
#else |
|
__mpfr_scale2 (d, exp) |
|
double d; |
|
int exp; |
|
#endif |
|
{ |
{ |
#if _GMP_IEEE_FLOATS |
int signd, sizetmp, inexact; |
{ |
unsigned int cnt, k; |
union ieee_double_extract x; |
mpfr_ptr tmp; |
x.d = d; |
TMP_DECL(marker); |
exp += x.s.exp; |
|
x.s.exp = exp; |
|
if (exp >= 2047) |
|
{ |
|
/* Return +-infinity */ |
|
x.s.exp = 2047; |
|
x.s.manl = x.s.manh = 0; |
|
} |
|
else if (exp < 1) |
|
{ |
|
x.s.exp = 1; /* smallest exponent (biased) */ |
|
/* Divide result by 2 until we have scaled it to the right IEEE |
|
denormalized number, but stop if it becomes zero. */ |
|
while (exp < 1 && x.d != 0) |
|
{ |
|
x.d *= 0.5; |
|
exp++; |
|
} |
|
} |
|
return x.d; |
|
} |
|
#else |
|
{ |
|
double factor, r; |
|
|
|
factor = 2.0; |
TMP_MARK(marker); |
if (exp < 0) |
MPFR_CLEAR_FLAGS(r); |
{ |
|
factor = 0.5; |
|
exp = -exp; |
|
} |
|
r = d; |
|
if (exp != 0) |
|
{ |
|
if ((exp & 1) != 0) |
|
r *= factor; |
|
exp >>= 1; |
|
while (exp != 0) |
|
{ |
|
factor *= factor; |
|
if ((exp & 1) != 0) |
|
r *= factor; |
|
exp >>= 1; |
|
} |
|
} |
|
return r; |
|
} |
|
#endif |
|
} |
|
|
|
|
if (d == 0) |
|
{ |
|
union ieee_double_extract x; |
|
|
/* End of part included from gmp */ |
MPFR_SET_ZERO(r); |
|
/* set correct sign */ |
|
x.d = d; |
|
if (((x.s.sig == 1) && (MPFR_SIGN(r) > 0)) |
|
|| ((x.s.sig == 0) && (MPFR_SIGN(r) < 0))) |
|
MPFR_CHANGE_SIGN(r); |
|
return 0; /* 0 is exact */ |
|
} |
|
else if (DOUBLE_ISNAN(d)) |
|
{ |
|
MPFR_SET_NAN(r); |
|
MPFR_RET_NAN; |
|
} |
|
else if (DOUBLE_ISINF(d)) |
|
{ |
|
MPFR_SET_INF(r); |
|
if ((d > 0 && (MPFR_SIGN(r) == -1)) || (d < 0 && (MPFR_SIGN(r) == 1))) |
|
MPFR_CHANGE_SIGN(r); |
|
return 0; /* infinity is exact */ |
|
} |
|
|
void |
/* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE, |
#if __STDC__ |
since PREC(r) may be different from PREC(tmp), and then both variables |
mpfr_set_d(mpfr_t r, double d, unsigned char rnd_mode) |
would have same precision in the mpfr_set4 call below. */ |
#else |
tmp = (mpfr_ptr) TMP_ALLOC(sizeof(mpfr_t)); |
mpfr_set_d(r, d, rnd_mode) |
MPFR_MANT(tmp) = TMP_ALLOC(MPFR_LIMBS_PER_DOUBLE * BYTES_PER_MP_LIMB); |
mpfr_t r; |
MPFR_PREC(tmp) = 53; |
double d; |
MPFR_SIZE(tmp) = MPFR_LIMBS_PER_DOUBLE; |
unsigned char rnd_mode; |
sizetmp = MPFR_LIMBS_PER_DOUBLE; |
#endif |
|
{ |
|
int signd, sizer; unsigned int cnt; |
|
|
|
if (d == 0) { SET_ZERO(r); return; } |
|
else if (isnan(d)) { SET_NAN(r); return; } |
|
|
|
signd = (d < 0) ? -1 : 1; |
signd = (d < 0) ? -1 : 1; |
d = ABS (d); |
d = ABS (d); |
sizer = (PREC(r)-1)/BITS_PER_MP_LIMB + 1; |
|
|
|
/* warning: __mpfr_extract_double requires at least two limbs */ |
MPFR_EXP(tmp) = __mpfr_extract_double (MPFR_MANT(tmp), d); |
if (sizer < MPFR_LIMBS_PER_DOUBLE) |
|
EXP(r) = __mpfr_extract_double (MANT(r), d, 0); |
|
else |
|
EXP(r) = __mpfr_extract_double (MANT(r) + sizer - MPFR_LIMBS_PER_DOUBLE, d, 1); |
|
|
|
if (sizer > MPFR_LIMBS_PER_DOUBLE) |
|
MPN_ZERO(MANT(r), sizer - MPFR_LIMBS_PER_DOUBLE); |
|
|
|
count_leading_zeros(cnt, MANT(r)[sizer-1]); |
/* determine number k of zero high limbs */ |
if (cnt) mpn_lshift(MANT(r), MANT(r), sizer, cnt); |
for (k = 0; k < sizetmp && MPFR_MANT(tmp)[sizetmp - 1 - k] == 0; k++); |
|
|
EXP(r) -= cnt; |
|
if (SIZE(r)*signd<0) CHANGE_SIGN(r); |
|
|
|
mpfr_round(r, rnd_mode, PREC(r)); |
count_leading_zeros (cnt, MPFR_MANT(tmp)[sizetmp - 1 - k]); |
return; |
|
} |
|
|
|
double |
if (cnt) |
#if __STDC__ |
mpn_lshift (MPFR_MANT(tmp) + k, MPFR_MANT(tmp), sizetmp - k, cnt); |
mpfr_get_d2(mpfr_srcptr src, long e) |
else if (k) |
#else |
MPN_COPY (MPFR_MANT(tmp) + k, MPFR_MANT(tmp), sizetmp - k); |
mpfr_get_d2(src, e) |
if (k) |
mpfr_srcptr(src); |
MPN_ZERO (MPFR_MANT(tmp), k); |
long e; |
|
#endif |
|
{ |
|
double res; |
|
mp_size_t size, i, n_limbs_to_use; |
|
mp_ptr qp; |
|
int negative; |
|
|
|
if (FLAG_NAN(src)) { |
MPFR_EXP(tmp) -= cnt + k * BITS_PER_MP_LIMB; |
#ifdef DEBUG |
|
printf("recognized NaN\n"); |
|
#endif |
|
return NaN; } |
|
if (NOTZERO(src)==0) return 0.0; |
|
size = 1+(PREC(src)-1)/BITS_PER_MP_LIMB; |
|
qp = MANT(src); |
|
negative = (SIGN(src)==-1); |
|
|
|
/* Warning: don't compute the abs(res) and set the sign afterwards, |
/* tmp is exact since PREC(tmp)=53 */ |
otherwise the current machine rounding mode will not be taken |
inexact = mpfr_set4 (r, tmp, rnd_mode, signd); |
correctly into account. */ |
|
/* res = (negative) ? -(double)qp[size - 1] : qp[size - 1]; */ |
|
res = 0.0; |
|
/* Warning: an arbitrary number of limbs may be required for an exact |
|
rounding. The following code is correct but not optimal since one |
|
may be able to decide without considering all limbs. */ |
|
/* n_limbs_to_use = MIN (MPFR_LIMBS_PER_DOUBLE, size); */ |
|
n_limbs_to_use = size; |
|
/* Accumulate the limbs from less significant to most significant |
|
otherwise due to rounding we may accumulate several ulps, |
|
especially in rounding towards -/+infinity. */ |
|
for (i = n_limbs_to_use; i>=1; i--) |
|
res = res / MP_BASE_AS_DOUBLE + |
|
((negative) ? -(double)qp[size - i] : qp[size - i]); |
|
res = __mpfr_scale2 (res, e - BITS_PER_MP_LIMB); |
|
|
|
return res; |
TMP_FREE(marker); |
|
return inexact; |
} |
} |
|
|
double |
|
#if __STDC__ |
|
mpfr_get_d(mpfr_srcptr src) |
|
#else |
|
mpfr_get_d(src) |
|
mpfr_srcptr src; |
|
#endif |
|
{ |
|
return mpfr_get_d2(src, EXP(src)); |
|
} |
|
|
|