1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
/* mpfr_get_ld -- convert a multiple precision floating-point number
to a machine long double
Copyright 2002, 2003 Free Software Foundation, Inc.
This file is part of the MPFR Library.
The MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
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
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include <stdio.h>
#include <float.h>
#include "gmp.h"
#include "gmp-impl.h"
#include "mpfr.h"
#include "mpfr-impl.h"
#ifndef LDBL_MANT_DIG
#define LDBL_MANT_DIG 113 /* works also if long double == quad */
#endif
#ifndef DBL_MANT_DIG
#define DBL_MANT_DIG 53
#endif
/* Various i386 systems have been seen with incorrect LDBL constants in
float.h (notes in set_ld.c), so force the value we know is right for IEEE
extended. */
#if HAVE_LDOUBLE_IEEE_EXT_LITTLE
#define MPFR_LDBL_MANT_DIG 64
#else
#define MPFR_LDBL_MANT_DIG LDBL_MANT_DIG
#endif
long double
mpfr_get_ld (mpfr_srcptr x, mp_rnd_t rnd_mode)
{
if (!mpfr_number_p (x)) /* NaN or Inf: check before 0.0 */
{
return (long double) mpfr_get_d (x, rnd_mode);
}
else if (MPFR_IS_ZERO(x))
{
return MPFR_SIGN(x) < 0 ? -0.0 : 0.0;
}
else /* now x is a normal non-zero number */
{
long double r; /* result */
long double m;
double s; /* part of result */
mp_exp_t e; /* exponent of x */
mp_exp_t sh; /* exponent shift, so that x/2^sh is in the double range */
mpfr_t y, z;
/* 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_set (y, x, rnd_mode);
e = MPFR_GET_EXP (y);
if (e > 1023)
{
sh = e - 1023;
MPFR_SET_EXP (y, 1023);
}
else if (e < -1021)
{
sh = e + 1021;
MPFR_SET_EXP (y, -1021);
}
else
{
sh = 0;
}
/* now -1021 <= e - sh = EXP(y) <= 1023 */
r = 0.0;
mpfr_init2 (z, DBL_MANT_DIG);
do
{
s = mpfr_get_d (y, GMP_RNDN); /* high part of y */
r += (long double) s;
mpfr_set_d (z, s, rnd_mode); /* exact */
mpfr_sub (y, y, z, rnd_mode); /* exact */
}
while (MPFR_NOTZERO(y));
/* we now have to multiply back by 2^sh */
if (sh != 0)
{
if (sh > 0)
m = 2.0;
else
{
m = 0.5;
sh = -sh;
}
e = 1; /* invariant: m = 2^e */
for (;;)
{
if (sh % 2)
r = r * m;
sh >>= 1;
if (sh == 0)
break;
m = m * m;
e = e + e;
}
}
mpfr_clear (z);
mpfr_clear (y);
return r;
}
}
|
