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