xref: /netbsd-src/external/lgpl3/mpfr/dist/src/set.c (revision ba125506a622fe649968631a56eba5d42ff57863)
1 /* mpfr_set -- copy of a floating-point number
2 
3 Copyright 1999, 2001-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 #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   /* Contrary to other mpfr_set4 based functions (mpfr_abs, mpfr_neg, etc.),
73      do not detect the case a == b as there is no interest to call mpfr_set
74      in this case, so that it is very unlikely that the user calls it
75      with a == b (this is the reverse of what is assumed for the other
76      functions). */
77   return mpfr_set4 (a, b, rnd_mode, MPFR_SIGN (b));
78 }
79 
80 /* Set a to |b| */
81 #undef mpfr_abs
82 int
mpfr_abs(mpfr_ptr a,mpfr_srcptr b,mpfr_rnd_t rnd_mode)83 mpfr_abs (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode)
84 {
85   if (MPFR_UNLIKELY (a != b))
86     return mpfr_set4 (a, b, rnd_mode, MPFR_SIGN_POS);
87   else
88     {
89       MPFR_SET_POS (a);
90       if (MPFR_UNLIKELY (MPFR_IS_NAN (b)))
91         MPFR_RET_NAN;
92       else
93         MPFR_RET (0);
94     }
95 }
96 
97 /* Round (u, inex) into s with rounding mode rnd_mode, where inex is
98    the ternary value associated with u, which has been obtained using
99    the *same* rounding mode rnd_mode.
100    Assumes PREC(u) = 2*PREC(s).
101    The generic algorithm is the following:
102    1. inex2 = mpfr_set (s, u, rnd_mode);
103    2. if (inex2 != 0) return inex2; else return inex;
104       except in the double-rounding case: in MPFR_RNDN, when u is in the
105       middle of two consecutive PREC(s)-bit numbers, if inex and inex2
106       are both > 0 (resp. both < 0), we correct s to nextbelow(s) (resp.
107       nextabove(s)), and return the opposite of inex.
108    Note: this function can be called with rnd_mode == MPFR_RNDF, in
109    which case, the rounding direction and the returned ternary value
110    are unspecified. */
111 int
mpfr_set_1_2(mpfr_ptr s,mpfr_srcptr u,mpfr_rnd_t rnd_mode,int inex)112 mpfr_set_1_2 (mpfr_ptr s, mpfr_srcptr u, mpfr_rnd_t rnd_mode, int inex)
113 {
114   mpfr_prec_t p = MPFR_PREC(s);
115   mpfr_prec_t sh = GMP_NUMB_BITS - p;
116   mp_limb_t rb, sb;
117   mp_limb_t *sp = MPFR_MANT(s);
118   mp_limb_t *up = MPFR_MANT(u);
119   mp_limb_t mask;
120   int inex2;
121 
122   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(u)))
123     {
124       mpfr_set (s, u, rnd_mode);
125       return inex;
126     }
127 
128   MPFR_ASSERTD(MPFR_PREC(u) == 2 * p);
129 
130   if (p < GMP_NUMB_BITS)
131     {
132       mask = MPFR_LIMB_MASK(sh);
133 
134       if (MPFR_PREC(u) <= GMP_NUMB_BITS)
135         {
136           mp_limb_t u0 = up[0];
137 
138           /* it suffices to round (u0, inex) */
139           rb = u0 & (MPFR_LIMB_ONE << (sh - 1));
140           sb = (u0 & mask) ^ rb;
141           sp[0] = u0 & ~mask;
142         }
143       else
144         {
145           mp_limb_t u1 = up[1];
146 
147           /* we need to round (u1, u0, inex) */
148           mask = MPFR_LIMB_MASK(sh);
149           rb = u1 & (MPFR_LIMB_ONE << (sh - 1));
150           sb = ((u1 & mask) ^ rb) | up[0];
151           sp[0] = u1 & ~mask;
152         }
153 
154       inex2 = inex * MPFR_SIGN(u);
155       MPFR_SIGN(s) = MPFR_SIGN(u);
156       MPFR_EXP(s) = MPFR_EXP(u);
157 
158       /* in case inex2 > 0, the value of u is rounded away,
159          thus we need to subtract something from (u0, rb, sb):
160          (a) if sb is not zero, since the subtracted value is < 1, we can leave
161          sb as it is;
162          (b) if rb <> 0 and sb = 0: change to rb = 0 and sb = 1
163          (c) if rb = sb = 0: change to rb = 1 and sb = 1, and subtract 1 */
164       if (inex2 > 0)
165         {
166           if (rb && sb == 0)
167             {
168               rb = 0;
169               sb = 1;
170             }
171         }
172       else /* inex2 <= 0 */
173         sb |= inex;
174 
175       /* now rb, sb are the round and sticky bits, together with the value of
176          sp[0], except possibly in the case rb = sb = 0 and inex2 > 0 */
177       if (rb == 0 && sb == 0)
178         {
179           if (inex2 <= 0)
180             MPFR_RET(0);
181           else /* inex2 > 0 can only occur for RNDN and RNDA:
182                   RNDN: return sp[0] and inex
183                   RNDA: return sp[0] and inex */
184             MPFR_RET(inex);
185         }
186       else if (rnd_mode == MPFR_RNDN)
187         {
188           if (rb == 0 || (sb == 0 && (sp[0] & (MPFR_LIMB_ONE << sh)) == 0))
189             goto truncate;
190           else
191             goto add_one_ulp;
192         }
193       else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(s)))
194         {
195         truncate:
196           MPFR_RET(-MPFR_SIGN(s));
197         }
198       else /* round away from zero */
199         {
200         add_one_ulp:
201           sp[0] += MPFR_LIMB_ONE << sh;
202           if (MPFR_UNLIKELY(sp[0] == 0))
203             {
204               sp[0] = MPFR_LIMB_HIGHBIT;
205               if (MPFR_EXP(s) + 1 <= __gmpfr_emax)
206                 MPFR_SET_EXP (s, MPFR_EXP(s) + 1);
207               else /* overflow */
208                 return mpfr_overflow (s, rnd_mode, MPFR_SIGN(s));
209             }
210           MPFR_RET(MPFR_SIGN(s));
211         }
212     }
213 
214   /* general case PREC(s) >= GMP_NUMB_BITS */
215   inex2 = mpfr_set (s, u, rnd_mode);
216   /* Check the double-rounding case, i.e. with u = middle of two
217      consecutive PREC(s)-bit numbers, which is equivalent to u being
218      exactly representable on PREC(s) + 1 bits but not on PREC(s) bits.
219      Moreover, since PREC(u) = 2*PREC(s), u and s cannot be identical
220      (as pointers), thus u was not changed. */
221   if (rnd_mode == MPFR_RNDN && inex * inex2 > 0 &&
222       mpfr_min_prec (u) == p + 1)
223     {
224       if (inex > 0)
225         mpfr_nextbelow (s);
226       else
227         mpfr_nextabove (s);
228       return -inex;
229     }
230   return inex2 != 0 ? inex2 : inex;
231 }
232