w32tex
About: TeX Live provides a comprehensive TeX system including all the major TeX-related programs, macro packages, and fonts that are free software. Windows sources.
  Fossies Dox: w32tex-src.tar.xz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

rndna.c
Go to the documentation of this file.
1 /* mpfr_round_nearest_away -- round to nearest away
2 
3 Copyright 2012-2020 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5 
6 This file is part of the GNU MPFR Library.
7 
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17 
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22 
23 #include "mpfr-impl.h"
24 
25 /* Note: this doesn't work for 2^(emin-2). Currently, the best that can be
26  done is to extend the exponent range internally, and do not support the
27  case emin = MPFR_EMIN_MIN from the caller. */
28 
29 /* In order to work, we'll save the context within the mantissa of 'rop'.
30  For that, we need to do some low level stuff like for init2/clear to create
31  a mantissa of bigger size, which contains the extra context.
32  To add a new field of type T in the context, add its type in
33  mpfr_size_limb_extended_t (if it is not already present)
34  and add a new value for the enum mpfr_index_extended_t before MANTISSA. */
35 typedef union {
44 typedef enum {
53  MANTISSA
55 
56 #define MPFR_MALLOC_EXTENDED_SIZE(s) \
57  (MANTISSA * sizeof(mpfr_size_limb_extended_t) + \
58  MPFR_BYTES_PER_MP_LIMB * (size_t) (s))
59 
60 /* This function is called before the applied function
61  and prepares rop to give it one more bit of precision
62  and to save its old value within it. */
63 void
65 {
66  mpfr_t tmp;
69  mpfr_prec_t p;
70  int inexact;
71  MPFR_SAVE_EXPO_DECL (expo);
72 
73  /* when emin is the smallest possible value, we cannot
74  determine the correct round-to-nearest-away rounding for
75  0.25*2^emin, which gets rounded to 0 with nearest-even,
76  instead of 0.5^2^emin. */
78 
79  p = MPFR_PREC (rop) + 1;
80  MPFR_SAVE_EXPO_MARK (expo);
81 
82  /* We can't use mpfr_init2 since we need to store an additional context
83  within the mantissa. */
85  /* Allocate the context within the needed mantissa. */
89 
90  /* Save the context first. */
91  ext[ALLOC_SIZE].si = xsize;
92  ext[OLD_MANTISSA].pi = MPFR_MANT(rop);
93  ext[OLD_EXPONENT].ex = MPFR_EXP(rop);
94  ext[OLD_SIGN].sg = MPFR_SIGN(rop);
95  ext[OLD_PREC].pr = MPFR_PREC(rop);
96  ext[OLD_FLAGS].fl = expo.saved_flags;
97  ext[OLD_EXP_MIN].ex = expo.saved_emin;
98  ext[OLD_EXP_MAX].ex = expo.saved_emax;
99 
100  /* Create tmp as a proper NAN. */
101  MPFR_PREC(tmp) = p; /* Set prec */
102  MPFR_SET_POS(tmp); /* Set a sign */
103  MPFR_MANT(tmp) = (mp_limb_t*)(ext+MANTISSA); /* Set Mantissa ptr */
104  MPFR_SET_NAN(tmp); /* initializes to NaN */
105 
106  /* Note: This full initialization to NaN may be unnecessary because of
107  the mpfr_set just below, but this is cleaner in case internals would
108  change in the future (for instance, some fields of tmp could be read
109  before being set, and reading an uninitialized value is UB here). */
110 
111  /* Copy rop into tmp now (rop may be used as input later). */
112  MPFR_DBGRES (inexact = mpfr_set(tmp, rop, MPFR_RNDN));
113  MPFR_ASSERTD(inexact == 0); /* Shall be exact since prec(tmp) > prec(rop) */
114 
115  /* Overwrite rop with tmp. */
116  rop[0] = tmp[0];
117 
118  /* The new rop now contains a copy of the old rop, with one extra bit of
119  precision while keeping the old rop "hidden" within its mantissa. */
120 
121  return;
122 }
123 
124 /* This function is called after the applied function.
125  rop is the one prepared in the function above,
126  and contains the result of the applied function.
127  This function restores rop like it was before the
128  calls to mpfr_round_nearest_away_begin while
129  copying it back the result of the applied function
130  and performing additional roundings. */
131 int
133 {
134  mpfr_t tmp;
137  mpfr_prec_t n;
138  MPFR_SAVE_EXPO_DECL (expo);
139 
140  /* Get back the hidden context.
141  Note: The cast to void * prevents the compiler from emitting a warning
142  (or error), such as:
143  cast increases required alignment of target type
144  with the -Wcast-align GCC option. Correct alignment is a consequence
145  of the code that built rop.
146  */
147  ext = ((mpfr_size_limb_extended_t *) (void *) MPFR_MANT(rop)) - MANTISSA;
148 
149  /* Create tmp with the result of the function. */
150  tmp[0] = rop[0];
151 
152  /* Revert rop back to what it was before calling
153  mpfr_round_neareset_away_begin. */
154  MPFR_PREC(rop) = ext[OLD_PREC].pr;
155  MPFR_SIGN(rop) = ext[OLD_SIGN].sg;
156  MPFR_EXP(rop) = ext[OLD_EXPONENT].ex;
157  MPFR_MANT(rop) = ext[OLD_MANTISSA].pi;
158  MPFR_ASSERTD(MPFR_PREC(tmp) == MPFR_PREC(rop)+1);
159 
160  /* Restore the saved context. */
161  expo.saved_flags = ext[OLD_FLAGS].fl;
162  expo.saved_emin = ext[OLD_EXP_MIN].ex;
163  expo.saved_emax = ext[OLD_EXP_MAX].ex;
164  xsize = ext[ALLOC_SIZE].si;
165 
166  /* Perform RNDNA. */
167  n = MPFR_PREC(rop);
168  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp)))
169  mpfr_set (rop, tmp, MPFR_RNDN); /* inex unchanged */
170  else
171  {
172  int lastbit, sh;
173 
174  MPFR_UNSIGNED_MINUS_MODULO(sh, n + 1);
175  lastbit = (MPFR_MANT(tmp)[0] >> sh) & 1;
176 
177  if (lastbit == 0)
178  mpfr_set (rop, tmp, MPFR_RNDN); /* exact, inex unchanged */
179  else if (inex == 0) /* midpoint: round away from zero */
180  inex = mpfr_set (rop, tmp, MPFR_RNDA);
181  else /* lastbit == 1, inex != 0: double rounding */
182  inex = mpfr_set (rop, tmp, (inex > 0) ? MPFR_RNDD : MPFR_RNDU);
183  }
184 
186  MPFR_SAVE_EXPO_FREE (expo);
187 
188  /* special treatment for the case rop = +/- 2^(emin-2), which should be
189  rounded to +/- 2^(emin-1). We do as if it was rounded to zero, thus the
190  mpfr_check_range() call will round it to +/- 2^(emin-1). */
191  if (inex == 0 && mpfr_cmp_si_2exp (rop, (mpfr_sgn (rop) > 0) ? 1 : -1,
192  __gmpfr_emin - 2) == 0)
193  inex = -mpfr_sgn (rop);
194 
195  /* Free tmp (cannot call mpfr_clear): free the associated context. */
197 
198  return mpfr_check_range (rop, inex, MPFR_RNDN);
199 }
Definition: asl.h:63
#define n
Definition: t4ht.c:1290
mpfr_exp_t __gmpfr_emin
Definition: exceptions.c:26
mpfr_flags_t __gmpfr_flags
Definition: exceptions.c:25
long int mp_size_t
Definition: gmp.h:175
small capitals from c petite p
Definition: afcover.h:72
void mpfr_free_func(void *ptr, size_t size)
Definition: mpfr-gmp.c:330
void * mpfr_allocate_func(size_t alloc_size)
Definition: mpfr-gmp.c:308
#define MPFR_SET_NAN(x)
Definition: mpfr-impl.h:1081
#define MPFR_UNLIKELY(x)
Definition: mpfr-impl.h:1490
#define MPFR_EMIN_MIN
Definition: mpfr-impl.h:1031
#define MPFR_SAVE_EXPO_UPDATE_FLAGS(x, flags)
Definition: mpfr-impl.h:1740
#define MPFR_SAVE_EXPO_MARK(x)
Definition: mpfr-impl.h:1730
#define MPFR_PREC(x)
Definition: mpfr-impl.h:958
#define MPFR_SET_POS(x)
Definition: mpfr-impl.h:1143
#define MPFR_DBGRES(A)
Definition: mpfr-impl.h:517
#define MPFR_SAVE_EXPO_FREE(x)
Definition: mpfr-impl.h:1736
#define mpfr_check_range(x, t, r)
Definition: mpfr-impl.h:1744
#define MPFR_SAVE_EXPO_DECL(x)
Definition: mpfr-impl.h:1729
#define MPFR_ASSERTN(expr)
Definition: mpfr-impl.h:495
#define MPFR_PREC2LIMBS(p)
Definition: mpfr-impl.h:955
#define MPFR_ASSERTD(expr)
Definition: mpfr-impl.h:516
#define MPFR_UNSIGNED_MINUS_MODULO(s, a)
Definition: mpfr-impl.h:1603
#define MPFR_EXP(x)
Definition: mpfr-impl.h:959
#define MPFR_IS_SINGULAR(x)
Definition: mpfr-impl.h:1100
#define MPFR_MANT(x)
Definition: mpfr-impl.h:960
int mpfr_cmp_si_2exp(mpfr_srcptr b, long int i, mpfr_exp_t f)
Definition: cmp_si.c:34
long mpfr_prec_t
Definition: mpfr.h:168
#define mpfr_set(a, b, r)
Definition: mpfr.h:864
__mpfr_struct mpfr_t[1]
Definition: mpfr.h:250
#define MPFR_PREC_MAX
Definition: mpfr.h:181
#define mpfr_sgn(_x)
Definition: mpfr.h:841
#define MPFR_SIGN(x)
Definition: mpfr.h:260
unsigned int mpfr_flags_t
Definition: mpfr.h:74
@ MPFR_RNDA
Definition: mpfr.h:107
@ MPFR_RNDN
Definition: mpfr.h:103
@ MPFR_RNDU
Definition: mpfr.h:105
@ MPFR_RNDD
Definition: mpfr.h:106
long mpfr_exp_t
Definition: mpfr.h:195
int mpfr_sign_t
Definition: mpfr.h:184
static int xsize
Definition: ppmtopjxl.c:41
void mpfr_round_nearest_away_begin(mpfr_t rop)
Definition: rndna.c:64
mpfr_index_extended_t
Definition: rndna.c:44
@ OLD_FLAGS
Definition: rndna.c:50
@ OLD_EXP_MIN
Definition: rndna.c:51
@ OLD_PREC
Definition: rndna.c:49
@ OLD_SIGN
Definition: rndna.c:48
@ MANTISSA
Definition: rndna.c:53
@ OLD_MANTISSA
Definition: rndna.c:46
@ OLD_EXPONENT
Definition: rndna.c:47
@ ALLOC_SIZE
Definition: rndna.c:45
@ OLD_EXP_MAX
Definition: rndna.c:52
int mpfr_round_nearest_away_end(mpfr_t rop, int inex)
Definition: rndna.c:132
#define MPFR_MALLOC_EXTENDED_SIZE(s)
Definition: rndna.c:56
char * ext
Definition: t4ht.c:938
mpfr_prec_t pr
Definition: rndna.c:39
mpfr_flags_t fl
Definition: rndna.c:41
mpfr_sign_t sg
Definition: rndna.c:40
mpfr_limb_ptr pi
Definition: rndna.c:42