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 326