xref: /netbsd-src/external/lgpl3/mpfr/dist/src/add1sp1_extracted.c (revision ba125506a622fe649968631a56eba5d42ff57863)
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
MPFR_Add1sp1_mk_state(int64_t sh,int64_t bx,uint64_t rb,uint64_t sb)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
mpfr_add1sp1(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode,mpfr_prec_t p)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 }
364