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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
|
/* mpfr_round_nearest_away -- round to nearest away
Copyright 2012-2017 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
The GNU 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 3 of the License, or (at your
option) any later version.
The GNU 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 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 "mpfr-impl.h"
/* Note: this doesn't work for 2^(emin-2). Currently, the best that can be
done is to extend the exponent range internally, and do not support the
case emin = MPFR_EMIN_MIN from the caller. */
/* In order to work, we'll save the context within the mantissa of 'rop'.
For that, we need to do some low level stuff like for init2/clear to create
a mantissa of bigger size, which contains the extra context.
To add a new field of type T in the context, add its type in
mpfr_size_limb_extended_t (if it is not already present)
and add a new value for the enum mpfr_index_extended_t before MANTISSA. */
typedef union {
mp_size_t si;
mp_limb_t li;
mpfr_exp_t ex;
mpfr_prec_t pr;
mpfr_sign_t sg;
mpfr_flags_t fl;
mpfr_limb_ptr pi;
} mpfr_size_limb_extended_t;
typedef enum {
ALLOC_SIZE = 0,
OLD_MANTISSA,
OLD_EXPONENT,
OLD_SIGN,
OLD_PREC,
OLD_FLAGS,
OLD_EXP_MIN,
OLD_EXP_MAX,
MANTISSA
} mpfr_index_extended_t ;
#define MPFR_MALLOC_EXTENDED_SIZE(s) \
(MANTISSA * sizeof(mpfr_size_limb_extended_t) + \
MPFR_BYTES_PER_MP_LIMB * (size_t) (s))
/* This function is called before the applied function
and prepares rop to give it one more bit of precision
and to save its old value within it. */
void
mpfr_round_nearest_away_begin (mpfr_t rop)
{
mpfr_t tmp;
mp_size_t xsize;
mpfr_size_limb_extended_t *ext;
mpfr_prec_t p;
int inexact;
MPFR_SAVE_EXPO_DECL (expo);
/* when emin is the smallest possible value, we cannot
determine the correct round-to-nearest-away rounding for
0.25*2^emin, which gets rounded to 0 with nearest-even,
instead of 0.5^2^emin. */
MPFR_ASSERTN(__gmpfr_emin > MPFR_EMIN_MIN);
p = MPFR_PREC (rop) + 1;
MPFR_SAVE_EXPO_MARK (expo);
/* We can't use mpfr_init2 since we need to store an additional context
within the mantissa. */
MPFR_ASSERTN(p <= MPFR_PREC_MAX);
/* Allocate the context within the needed mantissa. */
xsize = MPFR_PREC2LIMBS (p);
ext = (mpfr_size_limb_extended_t *)
mpfr_allocate_func (MPFR_MALLOC_EXTENDED_SIZE(xsize));
/* Save the context first. */
ext[ALLOC_SIZE].si = xsize;
ext[OLD_MANTISSA].pi = MPFR_MANT(rop);
ext[OLD_EXPONENT].ex = MPFR_EXP(rop);
ext[OLD_SIGN].sg = MPFR_SIGN(rop);
ext[OLD_PREC].pr = MPFR_PREC(rop);
ext[OLD_FLAGS].fl = expo.saved_flags;
ext[OLD_EXP_MIN].ex = expo.saved_emin;
ext[OLD_EXP_MAX].ex = expo.saved_emax;
/* Create tmp as a proper NAN. */
MPFR_PREC(tmp) = p; /* Set prec */
MPFR_SET_POS(tmp); /* Set a sign */
MPFR_MANT(tmp) = (mp_limb_t*)(ext+MANTISSA); /* Set Mantissa ptr */
MPFR_SET_NAN(tmp); /* initializes to NaN */
/* Note: This full initialization to NaN may be unnecessary because of
the mpfr_set just below, but this is cleaner in case internals would
change in the future (for instance, some fields of tmp could be read
before being set, and reading an uninitialized value is UB here). */
/* Copy rop into tmp now (rop may be used as input later). */
MPFR_DBGRES (inexact = mpfr_set(tmp, rop, MPFR_RNDN));
MPFR_ASSERTD(inexact == 0); /* Shall be exact since prec(tmp) > prec(rop) */
/* Overwrite rop with tmp. */
rop[0] = tmp[0];
/* The new rop now contains a copy of the old rop, with one extra bit of
precision while keeping the old rop "hidden" within its mantissa. */
return;
}
/* This function is called after the applied function.
rop is the one prepared in the function above,
and contains the result of the applied function.
This function restores rop like it was before the
calls to mpfr_round_nearest_away_begin while
copying it back the result of the applied function
and performing additional roundings. */
int
mpfr_round_nearest_away_end (mpfr_t rop, int inex)
{
mpfr_t tmp;
mp_size_t xsize;
mpfr_size_limb_extended_t *ext;
mpfr_prec_t n;
MPFR_SAVE_EXPO_DECL (expo);
/* Get back the hidden context. */
ext = ((mpfr_size_limb_extended_t *) MPFR_MANT(rop)) - MANTISSA;
/* Create tmp with the result of the function. */
tmp[0] = rop[0];
/* Revert rop back to what it was before calling
mpfr_round_neareset_away_begin. */
MPFR_PREC(rop) = ext[OLD_PREC].pr;
MPFR_SIGN(rop) = ext[OLD_SIGN].sg;
MPFR_EXP(rop) = ext[OLD_EXPONENT].ex;
MPFR_MANT(rop) = ext[OLD_MANTISSA].pi;
MPFR_ASSERTD(MPFR_PREC(tmp) == MPFR_PREC(rop)+1);
/* Restore the saved context. */
expo.saved_flags = ext[OLD_FLAGS].fl;
expo.saved_emin = ext[OLD_EXP_MIN].ex;
expo.saved_emax = ext[OLD_EXP_MAX].ex;
xsize = ext[ALLOC_SIZE].si;
/* Perform RNDNA. */
n = MPFR_PREC(rop);
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp)))
mpfr_set (rop, tmp, MPFR_RNDN); /* inex unchanged */
else
{
int lastbit, sh;
MPFR_UNSIGNED_MINUS_MODULO(sh, n + 1);
lastbit = (MPFR_MANT(tmp)[0] >> sh) & 1;
if (lastbit == 0)
mpfr_set (rop, tmp, MPFR_RNDN); /* exact, inex unchanged */
else if (inex == 0) /* midpoint: round away from zero */
inex = mpfr_set (rop, tmp, MPFR_RNDA);
else /* lastbit == 1, inex != 0: double rounding */
inex = mpfr_set (rop, tmp, (inex > 0) ? MPFR_RNDD : MPFR_RNDU);
}
MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
MPFR_SAVE_EXPO_FREE (expo);
/* special treatment for the case rop = +/- 2^(emin-2), which should be
rounded to +/- 2^(emin-1). We do as if it was rounded to zero, thus the
mpfr_check_range() call will round it to +/- 2^(emin-1). */
if (inex == 0 && mpfr_cmp_si_2exp (rop, (mpfr_sgn (rop) > 0) ? 1 : -1,
__gmpfr_emin - 2) == 0)
inex = -mpfr_sgn (rop);
/* Free tmp (cannot call mpfr_clear): free the associated context. */
mpfr_free_func(ext, MPFR_MALLOC_EXTENDED_SIZE(xsize));
return mpfr_check_range (rop, inex, MPFR_RNDN);
}
|
