1 /* mpfr_set -- copy of a floating-point number
2 
3 Copyright 1999, 2001-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 /* set a to abs(b) * signb: a=b when signb = SIGN(b), a=abs(b) when signb=1 */
26 MPFR_HOT_FUNCTION_ATTR int
mpfr_set4(mpfr_ptr a,mpfr_srcptr b,mpfr_rnd_t rnd_mode,int signb)27 mpfr_set4 (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode, int signb)
28 {
29   /* Sign is ALWAYS copied */
30   MPFR_SET_SIGN (a, signb);
31 
32   /* Exponent is also always copied since if the number is singular,
33      the exponent field determined the number.
34      Can't use MPFR_SET_EXP since the exponent may be singular */
35   MPFR_EXP (a) = MPFR_EXP (b);
36 
37   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (b)))
38     {
39       /* MPFR_SET_NAN, MPFR_SET_ZERO and MPFR_SET_INF are useless
40          since MPFR_EXP (a) = MPFR_EXP (b) does the job */
41       if (MPFR_IS_NAN (b))
42         MPFR_RET_NAN;
43       else
44         MPFR_RET (0);
45     }
46   else if (MPFR_PREC (b) == MPFR_PREC (a))
47     {
48       /* Same precision and b is not singular:
49        * just copy the mantissa, and set the exponent and the sign
50        * The result is exact. */
51       MPN_COPY (MPFR_MANT (a), MPFR_MANT (b), MPFR_LIMB_SIZE (b));
52       MPFR_RET (0);
53     }
54   else
55     {
56       int inex;
57 
58       /* Else Round B inside a */
59       MPFR_RNDRAW (inex, a, MPFR_MANT (b), MPFR_PREC (b), rnd_mode, signb,
60                    if (MPFR_UNLIKELY (++ MPFR_EXP (a) > __gmpfr_emax))
61                      return mpfr_overflow (a, rnd_mode, signb)
62                    );
63       MPFR_RET (inex);
64     }
65 }
66 
67 /* Set a to b  */
68 #undef mpfr_set
69 int
mpfr_set(mpfr_ptr a,mpfr_srcptr b,mpfr_rnd_t rnd_mode)70 mpfr_set (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
71 {
72   return mpfr_set4 (a, b, rnd_mode, MPFR_SIGN (b));
73 }
74 
75 /* Set a to |b| */
76 #undef mpfr_abs
77 int
mpfr_abs(mpfr_ptr a,mpfr_srcptr b,mpfr_rnd_t rnd_mode)78 mpfr_abs (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
79 {
80   return mpfr_set4 (a, b, rnd_mode, MPFR_SIGN_POS);
81 }
82 
83 /* Round (u, inex) into s with rounding mode rnd_mode, where inex is
84    the ternary value associated with u, which has been obtained using
85    the *same* rounding mode rnd_mode.
86    Assumes PREC(u) = 2*PREC(s).
87    The generic algorithm is the following:
88    1. inex2 = mpfr_set (s, u, rnd_mode);
89    2. if (inex2 != 0) return inex2; else return inex;
90       except in the double-rounding case: in MPFR_RNDN, when u is in the
91       middle of two consecutive PREC(s)-bit numbers, if inex and inex2
92       are both > 0 (resp. both < 0), we correct s to nextbelow(s) (resp.
93       nextabove(s)), and return the opposite of inex.
94    Note: this function can be called with rnd_mode == MPFR_RNDF, in
95    which case, the rounding direction and the returned ternary value
96    are unspecified. */
97 int
mpfr_set_1_2(mpfr_ptr s,mpfr_srcptr u,mpfr_rnd_t rnd_mode,int inex)98 mpfr_set_1_2 (mpfr_ptr s, mpfr_srcptr u, mpfr_rnd_t rnd_mode, int inex)
99 {
100   mpfr_prec_t p = MPFR_PREC(s);
101   mpfr_prec_t sh = GMP_NUMB_BITS - p;
102   mp_limb_t rb, sb;
103   mp_limb_t *sp = MPFR_MANT(s);
104   mp_limb_t *up = MPFR_MANT(u);
105   mp_limb_t mask;
106   int inex2;
107 
108   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(u)))
109     {
110       mpfr_set (s, u, rnd_mode);
111       return inex;
112     }
113 
114   MPFR_ASSERTD(MPFR_PREC(u) == 2 * p);
115 
116   if (p < GMP_NUMB_BITS)
117     {
118       mask = MPFR_LIMB_MASK(sh);
119 
120       if (MPFR_PREC(u) <= GMP_NUMB_BITS)
121         {
122           mp_limb_t u0 = up[0];
123 
124           /* it suffices to round (u0, inex) */
125           rb = u0 & (MPFR_LIMB_ONE << (sh - 1));
126           sb = (u0 & mask) ^ rb;
127           sp[0] = u0 & ~mask;
128         }
129       else
130         {
131           mp_limb_t u1 = up[1];
132 
133           /* we need to round (u1, u0, inex) */
134           mask = MPFR_LIMB_MASK(sh);
135           rb = u1 & (MPFR_LIMB_ONE << (sh - 1));
136           sb = ((u1 & mask) ^ rb) | up[0];
137           sp[0] = u1 & ~mask;
138         }
139 
140       inex2 = inex * MPFR_SIGN(u);
141       MPFR_SIGN(s) = MPFR_SIGN(u);
142       MPFR_EXP(s) = MPFR_EXP(u);
143 
144       /* in case inex2 > 0, the value of u is rounded away,
145          thus we need to subtract something from (u0, rb, sb):
146          (a) if sb is not zero, since the subtracted value is < 1, we can leave
147          sb as it is;
148          (b) if rb <> 0 and sb = 0: change to rb = 0 and sb = 1
149          (c) if rb = sb = 0: change to rb = 1 and sb = 1, and subtract 1 */
150       if (inex2 > 0)
151         {
152           if (rb && sb == 0)
153             {
154               rb = 0;
155               sb = 1;
156             }
157         }
158       else /* inex2 <= 0 */
159         sb |= inex;
160 
161       /* now rb, sb are the round and sticky bits, together with the value of
162          sp[0], except possibly in the case rb = sb = 0 and inex2 > 0 */
163       if (rb == 0 && sb == 0)
164         {
165           if (inex2 <= 0)
166             MPFR_RET(0);
167           else /* inex2 > 0 can only occur for RNDN and RNDA:
168                   RNDN: return sp[0] and inex
169                   RNDA: return sp[0] and inex */
170             MPFR_RET(inex);
171         }
172       else if (rnd_mode == MPFR_RNDN)
173         {
174           if (rb == 0 || (sb == 0 && (sp[0] & (MPFR_LIMB_ONE << sh)) == 0))
175             goto truncate;
176           else
177             goto add_one_ulp;
178         }
179       else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(s)))
180         {
181         truncate:
182           MPFR_RET(-MPFR_SIGN(s));
183         }
184       else /* round away from zero */
185         {
186         add_one_ulp:
187           sp[0] += MPFR_LIMB_ONE << sh;
188           if (MPFR_UNLIKELY(sp[0] == 0))
189             {
190               sp[0] = MPFR_LIMB_HIGHBIT;
191               if (MPFR_EXP(s) + 1 <= __gmpfr_emax)
192                 MPFR_SET_EXP (s, MPFR_EXP(s) + 1);
193               else /* overflow */
194                 return mpfr_overflow (s, rnd_mode, MPFR_SIGN(s));
195             }
196           MPFR_RET(MPFR_SIGN(s));
197         }
198     }
199 
200   /* general case PREC(s) >= GMP_NUMB_BITS */
201   inex2 = mpfr_set (s, u, rnd_mode);
202   /* Check the double-rounding case, i.e. with u = middle of two
203      consecutive PREC(s)-bit numbers, which is equivalent to u being
204      exactly representable on PREC(s) + 1 bits but not on PREC(s) bits.
205      Moreover, since PREC(u) = 2*PREC(s), u and s cannot be identical
206      (as pointers), thus u was not changed. */
207   if (rnd_mode == MPFR_RNDN && inex * inex2 > 0 &&
208       mpfr_min_prec (u) == p + 1)
209     {
210       if (inex > 0)
211         mpfr_nextbelow (s);
212       else
213         mpfr_nextabove (s);
214       return -inex;
215     }
216   return inex2 != 0 ? inex2 : inex;
217 }
218