1 /*
2 Copyright (C) 2014 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "flint/double_extras.h"
13 #include "arf.h"
14
15 /* most double: (2^53-1) * 2^971 */
16 /* least normal: 2^-1022 */
17 /* least subnormal: 2^-1074 */
18
19 static double
huge_double(arf_rnd_t rnd,int negative)20 huge_double(arf_rnd_t rnd, int negative)
21 {
22 double v;
23
24 if (rnd == ARF_RND_NEAR || rounds_up(rnd, negative))
25 v = D_INF;
26 else
27 v = ldexp(9007199254740991.0, 971);
28
29 return negative ? -v : v;
30 }
31
32 static double
tiny_double(arf_rnd_t rnd,int negative)33 tiny_double(arf_rnd_t rnd, int negative)
34 {
35 double v;
36
37 if (rnd == ARF_RND_NEAR || !rounds_up(rnd, negative))
38 v = 0.0;
39 else
40 v = ldexp(1.0, -1074);
41
42 return negative ? -v : v;
43 }
44
45 double
arf_get_d(const arf_t x,arf_rnd_t rnd)46 arf_get_d(const arf_t x, arf_rnd_t rnd)
47 {
48 if (arf_is_special(x))
49 {
50 if (arf_is_zero(x))
51 return 0.0;
52 else if (arf_is_pos_inf(x))
53 return D_INF;
54 else if (arf_is_neg_inf(x))
55 return -D_INF;
56 else
57 return D_NAN;
58 }
59 else
60 {
61 arf_t t;
62 mp_srcptr tp;
63 mp_size_t tn;
64 double v;
65
66 /* also catches bignum exponents */
67 if (ARF_EXP(x) > 1030 || ARF_EXP(x) < -1080)
68 {
69 if (fmpz_sgn(ARF_EXPREF(x)) > 0)
70 return huge_double(rnd, ARF_SGNBIT(x));
71 else
72 return tiny_double(rnd, ARF_SGNBIT(x));
73 }
74
75 /* allow mpfr to take care of corner cases for now */
76 if (ARF_EXP(x) > 1020 || ARF_EXP(x) <= -1020 || rnd == ARF_RND_NEAR)
77 {
78 mpfr_t xx;
79 ARF_GET_MPN_READONLY(tp, tn, x);
80
81 xx->_mpfr_d = (mp_ptr) tp;
82 xx->_mpfr_prec = tn * FLINT_BITS;
83 xx->_mpfr_sign = ARF_SGNBIT(x) ? -1 : 1;
84 xx->_mpfr_exp = ARF_EXP(x);
85
86 return mpfr_get_d(xx, rnd_to_mpfr(rnd));
87 }
88
89 arf_init(t);
90 arf_set_round(t, x, 53, rnd);
91 ARF_GET_MPN_READONLY(tp, tn, t);
92
93 if (tn == 1)
94 v = (double)(tp[0]);
95 else
96 v = (double)(tp[1]) + (double)(tp[0]) * ldexp(1,-32);
97
98 v = ldexp(v, ARF_EXP(t) - FLINT_BITS);
99
100 if (ARF_SGNBIT(t))
101 v = -v;
102
103 arf_clear(t);
104
105 return v;
106 }
107 }
108
109