1 /* mpfr_mul_1 -- internal function to perform a one-limb multiplication
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 #define bool int
31
32 typedef struct K___int64_t_uint64_t_uint64_t_s
33 {
34 int64_t fst;
35 uint64_t snd;
36 uint64_t thd;
37 }
38 K___int64_t_uint64_t_uint64_t;
39
40 #define MPFR_Lib_mpfr_struct __mpfr_struct
41 #define MPFR_Lib_mpfr_RET(I) (I) != 0 ? ((__gmpfr_flags |= MPFR_FLAGS_INEXACT), (I)) : 0
42 #define MPFR_Exceptions_mpfr_overflow mpfr_overflow
43 #define MPFR_Exceptions_mpfr_underflow mpfr_underflow
44 #define mpfr_prec _mpfr_prec
45 #define mpfr_exp _mpfr_exp
46 #define mpfr_d _mpfr_d
47 #define mpfr_sign _mpfr_sign
48 #define MPFR_Lib_gmp_NUMB_BITS GMP_NUMB_BITS
49 #define MPFR_Lib_mpfr_EMAX __gmpfr_emax
50 #define MPFR_Lib_mpfr_EMIN __gmpfr_emin
51 #define MPFR_RoundingMode_MPFR_RNDZ MPFR_RNDZ
52
53 #define MPFR_RoundingMode_uu___is_MPFR_RNDN(rnd_mode) (rnd_mode == MPFR_RNDN)
54 #define MPFR_RoundingMode_mpfr_IS_LIKE_RNDZ MPFR_IS_LIKE_RNDZ
55 #define MPFR_RoundingMode_mpfr_IS_LIKE_RNDA MPFR_IS_LIKE_RNDA
56
57 /* Special code for prec(a) < GMP_NUMB_BITS and
58 prec(b), prec(c) <= GMP_NUMB_BITS.
59 Note: this code was copied in sqr.c, function mpfr_sqr_1 (this saves a few cycles
60 with respect to have this function exported). As a consequence, any change here
61 should be reported in mpfr_sqr_1. */
62 static int
63 mpfr_mul_1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
64 mpfr_prec_t p)
65 {
66 uint64_t *ap = a->mpfr_d;
67 uint64_t *bp = b->mpfr_d;
68 uint64_t *cp = c->mpfr_d;
69 uint64_t b0 = bp[0U];
70 uint64_t c0 = cp[0U];
71 int64_t sh = MPFR_Lib_gmp_NUMB_BITS - p;
72 uint64_t mask = ((uint64_t)1U << (uint32_t)sh) - (uint64_t)1U;
73 int64_t ax = b->mpfr_exp + c->mpfr_exp;
74 //K___uint64_t_uint64_t scrut0 = MPFR_Umul_ppmm_umul_ppmm(b0, c0);
75 //uint64_t a0 = scrut0.fst;
76 //uint64_t sb = scrut0.snd;
77 uint64_t a0, sb;
78 umul_ppmm (a0, sb, b0, c0);
79 K___int64_t_uint64_t_uint64_t scrut;
80 if (a0 < (uint64_t)0x8000000000000000U)
81 scrut =
82 (
83 (K___int64_t_uint64_t_uint64_t){
84 .fst = ax - (int64_t)1,
85 .snd = a0 << (uint32_t)1U | sb >> (uint32_t)(MPFR_Lib_gmp_NUMB_BITS - (int64_t)1),
86 .thd = sb << (uint32_t)1U
87 }
88 );
89 else
90 scrut = ((K___int64_t_uint64_t_uint64_t){ .fst = ax, .snd = a0, .thd = sb });
91 int64_t ax1 = scrut.fst;
92 uint64_t a01 = scrut.snd;
93 uint64_t sb1 = scrut.thd;
94 uint64_t rb = a01 & (uint64_t)1U << (uint32_t)(sh - (int64_t)1);
95 uint64_t sb2 = sb1 | ((a01 & mask) ^ rb);
96 ap[0U] = a01 & ~mask;
97 MPFR_Lib_mpfr_struct uu___73_2756 = a[0U];
98 a[0U] =
99 (
100 (MPFR_Lib_mpfr_struct){
101 .mpfr_prec = uu___73_2756.mpfr_prec,
102 .mpfr_sign = b->mpfr_sign * c->mpfr_sign,
103 .mpfr_exp = uu___73_2756.mpfr_exp,
104 .mpfr_d = uu___73_2756.mpfr_d
105 }
106 );
107 uint64_t *ap1 = a->mpfr_d;
108 uint64_t a02 = ap1[0U];
109 if (ax1 > MPFR_Lib_mpfr_EMAX)
110 return MPFR_Exceptions_mpfr_overflow(a, rnd_mode, a->mpfr_sign);
111 else if (ax1 < MPFR_Lib_mpfr_EMIN)
112 {
113 bool aneg = a->mpfr_sign < (int32_t)0;
114 if
115 (
116 ax1
117 == MPFR_Lib_mpfr_EMIN - (int64_t)1
118 && a02 == ~mask
119 &&
120 ((MPFR_RoundingMode_uu___is_MPFR_RNDN(rnd_mode) && rb > (uint64_t)0U)
121 || ((rb | sb2) > (uint64_t)0U && MPFR_RoundingMode_mpfr_IS_LIKE_RNDA(rnd_mode, aneg)))
122 )
123 {
124 uint64_t *ap2 = a->mpfr_d;
125 uint64_t a03 = ap2[0U];
126 MPFR_Lib_mpfr_struct uu___72_2907 = a[0U];
127 a[0U] =
128 (
129 (MPFR_Lib_mpfr_struct){
130 .mpfr_prec = uu___72_2907.mpfr_prec,
131 .mpfr_sign = uu___72_2907.mpfr_sign,
132 .mpfr_exp = ax1,
133 .mpfr_d = uu___72_2907.mpfr_d
134 }
135 );
136 if (rb == (uint64_t)0U && sb2 == (uint64_t)0U)
137 return MPFR_Lib_mpfr_RET((int32_t)0);
138 else if (MPFR_RoundingMode_uu___is_MPFR_RNDN(rnd_mode))
139 if
140 (
141 rb
142 == (uint64_t)0U
143 || (sb2 == (uint64_t)0U && (a03 & (uint64_t)1U << (uint32_t)sh) == (uint64_t)0U)
144 )
145 {
146 int32_t ite;
147 if (a->mpfr_sign == (int32_t)1)
148 ite = (int32_t)-1;
149 else
150 ite = (int32_t)1;
151 return MPFR_Lib_mpfr_RET(ite);
152 }
153 else
154 {
155 uint64_t *ap3 = a->mpfr_d;
156 ap3[0U] = ap3[0U] + ((uint64_t)1U << (uint32_t)sh);
157 if (ap3[0U] == (uint64_t)0U)
158 {
159 ap3[0U] = (uint64_t)0x8000000000000000U;
160 if (ax1 + (int64_t)1 > MPFR_Lib_mpfr_EMAX)
161 return MPFR_Exceptions_mpfr_overflow(a, rnd_mode, a->mpfr_sign);
162 else
163 {
164 MPFR_Lib_mpfr_struct uu___72_3047 = a[0U];
165 a[0U] =
166 (
167 (MPFR_Lib_mpfr_struct){
168 .mpfr_prec = uu___72_3047.mpfr_prec,
169 .mpfr_sign = uu___72_3047.mpfr_sign,
170 .mpfr_exp = ax1 + (int64_t)1,
171 .mpfr_d = uu___72_3047.mpfr_d
172 }
173 );
174 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
175 }
176 }
177 else
178 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
179 }
180 else if (MPFR_RoundingMode_mpfr_IS_LIKE_RNDZ(rnd_mode, a->mpfr_sign < (int32_t)0))
181 {
182 int32_t ite;
183 if (a->mpfr_sign == (int32_t)1)
184 ite = (int32_t)-1;
185 else
186 ite = (int32_t)1;
187 return MPFR_Lib_mpfr_RET(ite);
188 }
189 else
190 {
191 uint64_t *ap3 = a->mpfr_d;
192 ap3[0U] = ap3[0U] + ((uint64_t)1U << (uint32_t)sh);
193 if (ap3[0U] == (uint64_t)0U)
194 {
195 ap3[0U] = (uint64_t)0x8000000000000000U;
196 if (ax1 + (int64_t)1 > MPFR_Lib_mpfr_EMAX)
197 return MPFR_Exceptions_mpfr_overflow(a, rnd_mode, a->mpfr_sign);
198 else
199 {
200 MPFR_Lib_mpfr_struct uu___72_3237 = a[0U];
201 a[0U] =
202 (
203 (MPFR_Lib_mpfr_struct){
204 .mpfr_prec = uu___72_3237.mpfr_prec,
205 .mpfr_sign = uu___72_3237.mpfr_sign,
206 .mpfr_exp = ax1 + (int64_t)1,
207 .mpfr_d = uu___72_3237.mpfr_d
208 }
209 );
210 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
211 }
212 }
213 else
214 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
215 }
216 }
217 else if
218 (
219 MPFR_RoundingMode_uu___is_MPFR_RNDN(rnd_mode)
220 &&
221 (ax1
222 < MPFR_Lib_mpfr_EMIN - (int64_t)1
223 || (a02 == (uint64_t)0x8000000000000000U && (rb | sb2) == (uint64_t)0U))
224 )
225 return MPFR_Exceptions_mpfr_underflow(a, MPFR_RoundingMode_MPFR_RNDZ, a->mpfr_sign);
226 else
227 return MPFR_Exceptions_mpfr_underflow(a, rnd_mode, a->mpfr_sign);
228 }
229 else
230 {
231 uint64_t *ap2 = a->mpfr_d;
232 uint64_t a03 = ap2[0U];
233 MPFR_Lib_mpfr_struct uu___72_3422 = a[0U];
234 a[0U] =
235 (
236 (MPFR_Lib_mpfr_struct){
237 .mpfr_prec = uu___72_3422.mpfr_prec,
238 .mpfr_sign = uu___72_3422.mpfr_sign,
239 .mpfr_exp = ax1,
240 .mpfr_d = uu___72_3422.mpfr_d
241 }
242 );
243 if (rb == (uint64_t)0U && sb2 == (uint64_t)0U)
244 return MPFR_Lib_mpfr_RET((int32_t)0);
245 else if (MPFR_RoundingMode_uu___is_MPFR_RNDN(rnd_mode))
246 if
247 (
248 rb
249 == (uint64_t)0U
250 || (sb2 == (uint64_t)0U && (a03 & (uint64_t)1U << (uint32_t)sh) == (uint64_t)0U)
251 )
252 {
253 int32_t ite;
254 if (a->mpfr_sign == (int32_t)1)
255 ite = (int32_t)-1;
256 else
257 ite = (int32_t)1;
258 return MPFR_Lib_mpfr_RET(ite);
259 }
260 else
261 {
262 uint64_t *ap3 = a->mpfr_d;
263 ap3[0U] = ap3[0U] + ((uint64_t)1U << (uint32_t)sh);
264 if (ap3[0U] == (uint64_t)0U)
265 {
266 ap3[0U] = (uint64_t)0x8000000000000000U;
267 if (ax1 + (int64_t)1 > MPFR_Lib_mpfr_EMAX)
268 return MPFR_Exceptions_mpfr_overflow(a, rnd_mode, a->mpfr_sign);
269 else
270 {
271 MPFR_Lib_mpfr_struct uu___72_3562 = a[0U];
272 a[0U] =
273 (
274 (MPFR_Lib_mpfr_struct){
275 .mpfr_prec = uu___72_3562.mpfr_prec,
276 .mpfr_sign = uu___72_3562.mpfr_sign,
277 .mpfr_exp = ax1 + (int64_t)1,
278 .mpfr_d = uu___72_3562.mpfr_d
279 }
280 );
281 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
282 }
283 }
284 else
285 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
286 }
287 else if (MPFR_RoundingMode_mpfr_IS_LIKE_RNDZ(rnd_mode, a->mpfr_sign < (int32_t)0))
288 {
289 int32_t ite;
290 if (a->mpfr_sign == (int32_t)1)
291 ite = (int32_t)-1;
292 else
293 ite = (int32_t)1;
294 return MPFR_Lib_mpfr_RET(ite);
295 }
296 else
297 {
298 uint64_t *ap3 = a->mpfr_d;
299 ap3[0U] = ap3[0U] + ((uint64_t)1U << (uint32_t)sh);
300 if (ap3[0U] == (uint64_t)0U)
301 {
302 ap3[0U] = (uint64_t)0x8000000000000000U;
303 if (ax1 + (int64_t)1 > MPFR_Lib_mpfr_EMAX)
304 return MPFR_Exceptions_mpfr_overflow(a, rnd_mode, a->mpfr_sign);
305 else
306 {
307 MPFR_Lib_mpfr_struct uu___72_3752 = a[0U];
308 a[0U] =
309 (
310 (MPFR_Lib_mpfr_struct){
311 .mpfr_prec = uu___72_3752.mpfr_prec,
312 .mpfr_sign = uu___72_3752.mpfr_sign,
313 .mpfr_exp = ax1 + (int64_t)1,
314 .mpfr_d = uu___72_3752.mpfr_d
315 }
316 );
317 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
318 }
319 }
320 else
321 return MPFR_Lib_mpfr_RET(a->mpfr_sign);
322 }
323 }
324 }
325