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