1 /* mpfr_add1sp1 -- internal function to perform a "real" addition on one limb
2 This code was extracted by Kremlin from a formal proof in F*
3 done by Jianyang Pan in April-August 2018: do not modify it!
4
5 Source: https://github.com/project-everest/hacl-star/tree/dev_mpfr/code/mpfr
6
7 Copyright 2004-2023 Free Software Foundation, Inc.
8 Contributed by the AriC and Caramba projects, INRIA.
9
10 This file is part of the GNU MPFR Library.
11
12 The GNU MPFR Library is free software; you can redistribute it and/or modify
13 it under the terms of the GNU Lesser General Public License as published by
14 the Free Software Foundation; either version 3 of the License, or (at your
15 option) any later version.
16
17 The GNU MPFR Library is distributed in the hope that it will be useful, but
18 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
19 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
20 License for more details.
21
22 You should have received a copy of the GNU Lesser General Public License
23 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
24 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
25 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
26
27 #define int64_t long
28 #define uint32_t unsigned int
29 #define uint64_t mp_limb_t
30
31 typedef struct MPFR_Add1sp1_state_s
32 {
33 int64_t sh;
34 int64_t bx;
35 uint64_t rb;
36 uint64_t sb;
37 }
38 MPFR_Add1sp1_state;
39
40 static MPFR_Add1sp1_state
41 MPFR_Add1sp1_mk_state(int64_t sh, int64_t bx, uint64_t rb, uint64_t sb)
42 {
43 return ((MPFR_Add1sp1_state){ .sh = sh, .bx = bx, .rb = rb, .sb = sb });
44 }
45
46 typedef struct K___uint64_t_int64_t_s
47 {
48 uint64_t fst;
49 int64_t snd;
50 }
51 K___uint64_t_int64_t;
52
53 typedef struct K___uint64_t_uint64_t_int64_t_s
54 {
55 uint64_t fst;
56 uint64_t snd;
57 int64_t thd;
58 }
59 K___uint64_t_uint64_t_int64_t;
60
61 #define MPFR_Lib_mpfr_struct __mpfr_struct
62 #define MPFR_Lib_mpfr_RET(I) ((I) != 0 ? ((__gmpfr_flags |= MPFR_FLAGS_INEXACT), (I)) : 0)
63 #define MPFR_Exceptions_mpfr_overflow mpfr_overflow
64 #define mpfr_prec _mpfr_prec
65 #define mpfr_exp _mpfr_exp
66 #define mpfr_d _mpfr_d
67 #define mpfr_sign _mpfr_sign
68 #define MPFR_Lib_gmp_NUMB_BITS GMP_NUMB_BITS
69 #define MPFR_Lib_mpfr_EMAX __gmpfr_emax
70
71 #define MPFR_RoundingMode_uu___is_MPFR_RNDN(rnd_mode) (rnd_mode == MPFR_RNDN)
72 #define MPFR_RoundingMode_mpfr_IS_LIKE_RNDZ MPFR_IS_LIKE_RNDZ
73
74 /* same as mpfr_add1sp, but for p < GMP_NUMB_BITS */
75 static int
76 mpfr_add1sp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
77 mpfr_prec_t p)
78 {
79 MPFR_Lib_mpfr_struct a0 = a[0U];
80 MPFR_Lib_mpfr_struct b0 = b[0U];
81 MPFR_Lib_mpfr_struct c0 = c[0U];
82 int64_t bx = b0.mpfr_exp;
83 int64_t cx = c0.mpfr_exp;
84 int64_t sh = MPFR_Lib_gmp_NUMB_BITS - p;
85 MPFR_Add1sp1_state st;
86 if (bx == cx)
87 {
88 uint64_t *ap = a0.mpfr_d;
89 uint64_t *bp = b0.mpfr_d;
90 uint64_t *cp = c0.mpfr_d;
91 uint64_t a01 = (bp[0U] >> (uint32_t)1U) + (cp[0U] >> (uint32_t)1U);
92 int64_t bx1 = b0.mpfr_exp + (int64_t)1;
93 uint64_t rb = a01 & (uint64_t)1U << (uint32_t)(sh - (int64_t)1);
94 ap[0U] = a01 ^ rb;
95 uint64_t sb = (uint64_t)0U;
96 st = MPFR_Add1sp1_mk_state(sh, bx1, rb, sb);
97 }
98 else
99 {
100 MPFR_Add1sp1_state ite0;
101 if (bx > cx)
102 {
103 int64_t bx1 = b0.mpfr_exp;
104 int64_t cx1 = c0.mpfr_exp;
105 int64_t d = bx1 - cx1;
106 uint64_t mask = ((uint64_t)1U << (uint32_t)sh) - (uint64_t)1U;
107 MPFR_Add1sp1_state ite1;
108 if (d < sh)
109 {
110 uint64_t *ap = a0.mpfr_d;
111 uint64_t *bp = b0.mpfr_d;
112 uint64_t *cp = c0.mpfr_d;
113 int64_t bx2 = b0.mpfr_exp;
114 uint64_t a01 = bp[0U] + (cp[0U] >> (uint32_t)d);
115 K___uint64_t_int64_t scrut;
116 if (a01 < bp[0U])
117 scrut =
118 (
119 (K___uint64_t_int64_t){
120 .fst = (uint64_t)0x8000000000000000U | a01 >> (uint32_t)1U,
121 .snd = bx2 + (int64_t)1
122 }
123 );
124 else
125 scrut = ((K___uint64_t_int64_t){ .fst = a01, .snd = bx2 });
126 uint64_t a02 = scrut.fst;
127 int64_t bx3 = scrut.snd;
128 uint64_t rb = a02 & (uint64_t)1U << (uint32_t)(sh - (int64_t)1);
129 uint64_t sb = (a02 & mask) ^ rb;
130 ap[0U] = a02 & ~mask;
131 ite1 = MPFR_Add1sp1_mk_state(sh, bx3, rb, sb);
132 }
133 else
134 {
135 MPFR_Add1sp1_state ite;
136 if (d < MPFR_Lib_gmp_NUMB_BITS)
137 {
138 uint64_t *ap = a0.mpfr_d;
139 uint64_t *bp = b0.mpfr_d;
140 uint64_t *cp = c0.mpfr_d;
141 int64_t bx2 = b0.mpfr_exp;
142 uint64_t sb = cp[0U] << (uint32_t)(MPFR_Lib_gmp_NUMB_BITS - d);
143 uint64_t a01 = bp[0U] + (cp[0U] >> (uint32_t)d);
144 K___uint64_t_uint64_t_int64_t scrut;
145 if (a01 < bp[0U])
146 scrut =
147 (
148 (K___uint64_t_uint64_t_int64_t){
149 .fst = sb | (a01 & (uint64_t)1U),
150 .snd = (uint64_t)0x8000000000000000U | a01 >> (uint32_t)1U,
151 .thd = bx2 + (int64_t)1
152 }
153 );
154 else
155 scrut = ((K___uint64_t_uint64_t_int64_t){ .fst = sb, .snd = a01, .thd = bx2 });
156 uint64_t sb1 = scrut.fst;
157 uint64_t a02 = scrut.snd;
158 int64_t bx3 = scrut.thd;
159 uint64_t rb = a02 & (uint64_t)1U << (uint32_t)(sh - (int64_t)1);
160 uint64_t sb2 = sb1 | ((a02 & mask) ^ rb);
161 ap[0U] = a02 & ~mask;
162 ite = MPFR_Add1sp1_mk_state(sh, bx3, rb, sb2);
163 }
164 else
165 {
166 uint64_t *ap = a0.mpfr_d;
167 uint64_t *bp = b0.mpfr_d;
168 int64_t bx2 = b0.mpfr_exp;
169 ap[0U] = bp[0U];
170 uint64_t rb = (uint64_t)0U;
171 uint64_t sb = (uint64_t)1U;
172 ite = MPFR_Add1sp1_mk_state(sh, bx2, rb, sb);
173 }
174 ite1 = ite;
175 }
176 ite0 = ite1;
177 }
178 else
179 {
180 int64_t bx1 = c0.mpfr_exp;
181 int64_t cx1 = b0.mpfr_exp;
182 int64_t d = bx1 - cx1;
183 uint64_t mask = ((uint64_t)1U << (uint32_t)sh) - (uint64_t)1U;
184 MPFR_Add1sp1_state ite1;
185 if (d < sh)
186 {
187 uint64_t *ap = a0.mpfr_d;
188 uint64_t *bp = c0.mpfr_d;
189 uint64_t *cp = b0.mpfr_d;
190 int64_t bx2 = c0.mpfr_exp;
191 uint64_t a01 = bp[0U] + (cp[0U] >> (uint32_t)d);
192 K___uint64_t_int64_t scrut;
193 if (a01 < bp[0U])
194 scrut =
195 (
196 (K___uint64_t_int64_t){
197 .fst = (uint64_t)0x8000000000000000U | a01 >> (uint32_t)1U,
198 .snd = bx2 + (int64_t)1
199 }
200 );
201 else
202 scrut = ((K___uint64_t_int64_t){ .fst = a01, .snd = bx2 });
203 uint64_t a02 = scrut.fst;
204 int64_t bx3 = scrut.snd;
205 uint64_t rb = a02 & (uint64_t)1U << (uint32_t)(sh - (int64_t)1);
206 uint64_t sb = (a02 & mask) ^ rb;
207 ap[0U] = a02 & ~mask;
208 ite1 = MPFR_Add1sp1_mk_state(sh, bx3, rb, sb);
209 }
210 else
211 {
212 MPFR_Add1sp1_state ite;
213 if (d < MPFR_Lib_gmp_NUMB_BITS)
214 {
215 uint64_t *ap = a0.mpfr_d;
216 uint64_t *bp = c0.mpfr_d;
217 uint64_t *cp = b0.mpfr_d;
218 int64_t bx2 = c0.mpfr_exp;
219 uint64_t sb = cp[0U] << (uint32_t)(MPFR_Lib_gmp_NUMB_BITS - d);
220 uint64_t a01 = bp[0U] + (cp[0U] >> (uint32_t)d);
221 K___uint64_t_uint64_t_int64_t scrut;
222 if (a01 < bp[0U])
223 scrut =
224 (
225 (K___uint64_t_uint64_t_int64_t){
226 .fst = sb | (a01 & (uint64_t)1U),
227 .snd = (uint64_t)0x8000000000000000U | a01 >> (uint32_t)1U,
228 .thd = bx2 + (int64_t)1
229 }
230 );
231 else
232 scrut = ((K___uint64_t_uint64_t_int64_t){ .fst = sb, .snd = a01, .thd = bx2 });
233 uint64_t sb1 = scrut.fst;
234 uint64_t a02 = scrut.snd;
235 int64_t bx3 = scrut.thd;
236 uint64_t rb = a02 & (uint64_t)1U << (uint32_t)(sh - (int64_t)1);
237 uint64_t sb2 = sb1 | ((a02 & mask) ^ rb);
238 ap[0U] = a02 & ~mask;
239 ite = MPFR_Add1sp1_mk_state(sh, bx3, rb, sb2);
240 }
241 else
242 {
243 uint64_t *ap = a0.mpfr_d;
244 uint64_t *bp = c0.mpfr_d;
245 int64_t bx2 = c0.mpfr_exp;
246 ap[0U] = bp[0U];
247 uint64_t rb = (uint64_t)0U;
248 uint64_t sb = (uint64_t)1U;
249 ite = MPFR_Add1sp1_mk_state(sh, bx2, rb, sb);
250 }
251 ite1 = ite;
252 }
253 ite0 = ite1;
254 }
255 st = ite0;
256 }
257 if (st.bx > MPFR_Lib_mpfr_EMAX)
258 {
259 int32_t t = MPFR_Exceptions_mpfr_overflow(a, rnd_mode, a->mpfr_sign);
260 return t;
261 }
262 else
263 {
264 uint64_t *ap = a->mpfr_d;
265 uint64_t a01 = ap[0U];
266 MPFR_Lib_mpfr_struct uu___72_3478 = a[0U];
267 a[0U] =
268 (
269 (MPFR_Lib_mpfr_struct){
270 .mpfr_prec = uu___72_3478.mpfr_prec,
271 .mpfr_sign = uu___72_3478.mpfr_sign,
272 .mpfr_exp = st.bx,
273 .mpfr_d = uu___72_3478.mpfr_d
274 }
275 );
276 if (st.rb == (uint64_t)0U && st.sb == (uint64_t)0U)
277 return MPFR_Lib_mpfr_RET((int32_t)0);
278 else if (MPFR_RoundingMode_uu___is_MPFR_RNDN(rnd_mode))
279 if
280 (
281 st.rb
282 == (uint64_t)0U
283 || (st.sb == (uint64_t)0U && (a01 & (uint64_t)1U << (uint32_t)st.sh) == (uint64_t)0U)
284 )
285 {
286 int32_t ite;
287 if (a->mpfr_sign == (int32_t)1)
288 ite = (int32_t)-1;
289 else
290 ite = (int32_t)1;
291 return MPFR_Lib_mpfr_RET(ite);
292 }
293 else
294 {
295 uint64_t *ap1 = a->mpfr_d;
296 ap1[0U] = ap1[0U] + ((uint64_t)1U << (uint32_t)st.sh);
297 if (ap1[0U] == (uint64_t)0U)
298 {
299 ap1[0U] = (uint64_t)0x8000000000000000U;
300 if (st.bx + (int64_t)1 <= MPFR_Lib_mpfr_EMAX)
301 {
302 MPFR_Lib_mpfr_struct uu___72_3574 = a[0U];
303 a[0U] =
304 (
305 (MPFR_Lib_mpfr_struct){
306 .mpfr_prec = uu___72_3574.mpfr_prec,
307 .mpfr_sign = uu___72_3574.mpfr_sign,
308 .mpfr_exp = st.bx + (int64_t)1,
309 .mpfr_d = uu___72_3574.mpfr_d
310 }
311 );
312 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
313 }
314 else
315 {
316 int32_t t = MPFR_Exceptions_mpfr_overflow(a, rnd_mode, a->mpfr_sign);
317 return MPFR_Lib_mpfr_RET(t);
318 }
319 }
320 else
321 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
322 }
323 else if (MPFR_RoundingMode_mpfr_IS_LIKE_RNDZ(rnd_mode, a->mpfr_sign < (int32_t)0))
324 {
325 int32_t ite;
326 if (a->mpfr_sign == (int32_t)1)
327 ite = (int32_t)-1;
328 else
329 ite = (int32_t)1;
330 return MPFR_Lib_mpfr_RET(ite);
331 }
332 else
333 {
334 uint64_t *ap1 = a->mpfr_d;
335 ap1[0U] = ap1[0U] + ((uint64_t)1U << (uint32_t)st.sh);
336 if (ap1[0U] == (uint64_t)0U)
337 {
338 ap1[0U] = (uint64_t)0x8000000000000000U;
339 if (st.bx + (int64_t)1 <= MPFR_Lib_mpfr_EMAX)
340 {
341 MPFR_Lib_mpfr_struct uu___72_3781 = a[0U];
342 a[0U] =
343 (
344 (MPFR_Lib_mpfr_struct){
345 .mpfr_prec = uu___72_3781.mpfr_prec,
346 .mpfr_sign = uu___72_3781.mpfr_sign,
347 .mpfr_exp = st.bx + (int64_t)1,
348 .mpfr_d = uu___72_3781.mpfr_d
349 }
350 );
351 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
352 }
353 else
354 {
355 int32_t t = MPFR_Exceptions_mpfr_overflow(a, rnd_mode, a->mpfr_sign);
356 return MPFR_Lib_mpfr_RET(t);
357 }
358 }
359 else
360 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
361 }
362 }
363 }