xref: /dragonfly/contrib/mpfr/src/get_f.c (revision ed36d35d)
1 /* mpfr_get_f -- convert a MPFR number to a GNU MPF number
2 
3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel 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 http://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 /* Since MPFR-3.0, return the usual inexact value.
26    The erange flag is set if an error occurred in the conversion
27    (y is NaN, +Inf, or -Inf that have no equivalent in mpf)
28 */
29 int
30 mpfr_get_f (mpf_ptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
31 {
32   int inex;
33   mp_size_t sx, sy;
34   mpfr_prec_t precx, precy;
35   mp_limb_t *xp;
36   int sh;
37 
38   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(y)))
39     {
40       if (MPFR_IS_ZERO(y))
41         {
42           mpf_set_ui (x, 0);
43           return 0;
44         }
45       else if (MPFR_IS_NAN (y))
46         {
47           MPFR_SET_ERANGE ();
48           return 0;
49         }
50       else /* y is plus infinity (resp. minus infinity), set x to the maximum
51               value (resp. the minimum value) in precision PREC(x) */
52         {
53           int i;
54           mp_limb_t *xp;
55 
56           MPFR_SET_ERANGE ();
57 
58           /* To this day, [mp_exp_t] and mp_size_t are #defined as the same
59              type */
60           EXP (x) = MP_SIZE_T_MAX;
61 
62           sx = PREC (x);
63           SIZ (x) = sx;
64           xp = PTR (x);
65           for (i = 0; i < sx; i++)
66             xp[i] = MP_LIMB_T_MAX;
67 
68           if (MPFR_IS_POS (y))
69             return -1;
70           else
71             {
72               mpf_neg (x, x);
73               return +1;
74             }
75         }
76     }
77 
78   sx = PREC(x); /* number of limbs of the mantissa of x */
79 
80   precy = MPFR_PREC(y);
81   precx = (mpfr_prec_t) sx * GMP_NUMB_BITS;
82   sy = MPFR_LIMB_SIZE (y);
83 
84   xp = PTR (x);
85 
86   /* since mpf numbers are represented in base 2^GMP_NUMB_BITS,
87      we loose -EXP(y) % GMP_NUMB_BITS bits in the most significant limb */
88   sh = MPFR_GET_EXP(y) % GMP_NUMB_BITS;
89   sh = sh <= 0 ? - sh : GMP_NUMB_BITS - sh;
90   MPFR_ASSERTD (sh >= 0);
91   if (precy + sh <= precx) /* we can copy directly */
92     {
93       mp_size_t ds;
94 
95       MPFR_ASSERTN (sx >= sy);
96       ds = sx - sy;
97 
98       if (sh != 0)
99         {
100           mp_limb_t out;
101           out = mpn_rshift (xp + ds, MPFR_MANT(y), sy, sh);
102           MPFR_ASSERTN (ds > 0 || out == 0);
103           if (ds > 0)
104             xp[--ds] = out;
105         }
106       else
107         MPN_COPY (xp + ds, MPFR_MANT (y), sy);
108       if (ds > 0)
109         MPN_ZERO (xp, ds);
110       EXP(x) = (MPFR_GET_EXP(y) + sh) / GMP_NUMB_BITS;
111       inex = 0;
112     }
113   else /* we have to round to precx - sh bits */
114     {
115       mpfr_t z;
116       mp_size_t sz;
117 
118       /* Recall that precx = (mpfr_prec_t) sx * GMP_NUMB_BITS, thus removing
119          sh bits (sh < GMP_NUMB_BITSS) won't reduce the number of limbs. */
120       mpfr_init2 (z, precx - sh);
121       sz = MPFR_LIMB_SIZE (z);
122       MPFR_ASSERTN (sx == sz);
123 
124       inex = mpfr_set (z, y, rnd_mode);
125       /* warning, sh may change due to rounding, but then z is a power of two,
126          thus we can safely ignore its last bit which is 0 */
127       sh = MPFR_GET_EXP(z) % GMP_NUMB_BITS;
128       sh = sh <= 0 ? - sh : GMP_NUMB_BITS - sh;
129       MPFR_ASSERTD (sh >= 0);
130       if (sh != 0)
131         {
132           mp_limb_t out;
133           out = mpn_rshift (xp, MPFR_MANT(z), sz, sh);
134           /* If sh hasn't changed, it is the number of the non-significant
135              bits in the lowest limb of z. Therefore out == 0. */
136           MPFR_ASSERTD (out == 0);  (void) out; /* avoid a warning */
137         }
138       else
139         MPN_COPY (xp, MPFR_MANT(z), sz);
140       EXP(x) = (MPFR_GET_EXP(z) + sh) / GMP_NUMB_BITS;
141       mpfr_clear (z);
142     }
143 
144   /* set size and sign */
145   SIZ(x) = (MPFR_FROM_SIGN_TO_INT(MPFR_SIGN(y)) < 0) ? -sx : sx;
146 
147   return inex;
148 }
149