1 /* Copyright (C) 2007-2020 Free Software Foundation, Inc. 2 3 This file is part of GCC. 4 5 GCC is free software; you can redistribute it and/or modify it under 6 the terms of the GNU General Public License as published by the Free 7 Software Foundation; either version 3, or (at your option) any later 8 version. 9 10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11 WARRANTY; without even the implied warranty of MERCHANTABILITY or 12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13 for more details. 14 15 Under Section 7 of GPL version 3, you are granted additional 16 permissions described in the GCC Runtime Library Exception, version 17 3.1, as published by the Free Software Foundation. 18 19 You should have received a copy of the GNU General Public License and 20 a copy of the GCC Runtime Library Exception along with this program; 21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22 <http://www.gnu.org/licenses/>. */ 23 24 #define BID_128RES 25 #include "bid_internal.h" 26 #include "bid_sqrt_macros.h" 27 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 28 #include <fenv.h> 29 30 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT 31 #endif 32 33 BID128_FUNCTION_ARG1 (bid128_sqrt, x) 34 35 UINT256 M256, C256, C4, C8; 36 UINT128 CX, CX1, CX2, A10, S2, T128, TP128, CS, CSM, res; 37 UINT64 sign_x, Carry; 38 SINT64 D; 39 int_float fx, f64; 40 int exponent_x, bin_expon_cx; 41 int digits, scale, exponent_q; 42 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 43 fexcept_t binaryflags = 0; 44 #endif 45 46 // unpack arguments, check for NaN or Infinity 47 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) { 48 res.w[1] = CX.w[1]; 49 res.w[0] = CX.w[0]; 50 // NaN ? 51 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 52 #ifdef SET_STATUS_FLAGS 53 if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN 54 __set_status_flags (pfpsf, INVALID_EXCEPTION); 55 #endif 56 res.w[1] = CX.w[1] & QUIET_MASK64; 57 BID_RETURN (res); 58 } 59 // x is Infinity? 60 if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) { 61 res.w[1] = CX.w[1]; 62 if (sign_x) { 63 // -Inf, return NaN 64 res.w[1] = 0x7c00000000000000ull; 65 #ifdef SET_STATUS_FLAGS 66 __set_status_flags (pfpsf, INVALID_EXCEPTION); 67 #endif 68 } 69 BID_RETURN (res); 70 } 71 // x is 0 otherwise 72 73 res.w[1] = 74 sign_x | 75 ((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1) << 49); 76 res.w[0] = 0; 77 BID_RETURN (res); 78 } 79 if (sign_x) { 80 res.w[1] = 0x7c00000000000000ull; 81 res.w[0] = 0; 82 #ifdef SET_STATUS_FLAGS 83 __set_status_flags (pfpsf, INVALID_EXCEPTION); 84 #endif 85 BID_RETURN (res); 86 } 87 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 88 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS); 89 #endif 90 // 2^64 91 f64.i = 0x5f800000; 92 93 // fx ~ CX 94 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0]; 95 bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f; 96 digits = estimate_decimal_digits[bin_expon_cx]; 97 98 A10 = CX; 99 if (exponent_x & 1) { 100 A10.w[1] = (CX.w[1] << 3) | (CX.w[0] >> 61); 101 A10.w[0] = CX.w[0] << 3; 102 CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63); 103 CX2.w[0] = CX.w[0] << 1; 104 __add_128_128 (A10, A10, CX2); 105 } 106 107 CS.w[0] = short_sqrt128 (A10); 108 CS.w[1] = 0; 109 // check for exact result 110 if (CS.w[0] * CS.w[0] == A10.w[0]) { 111 __mul_64x64_to_128_fast (S2, CS.w[0], CS.w[0]); 112 if (S2.w[1] == A10.w[1]) // && S2.w[0]==A10.w[0]) 113 { 114 get_BID128_very_fast (&res, 0, 115 (exponent_x + 116 DECIMAL_EXPONENT_BIAS_128) >> 1, CS); 117 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 118 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 119 #endif 120 BID_RETURN (res); 121 } 122 } 123 // get number of digits in CX 124 D = CX.w[1] - power10_index_binexp_128[bin_expon_cx].w[1]; 125 if (D > 0 126 || (!D && CX.w[0] >= power10_index_binexp_128[bin_expon_cx].w[0])) 127 digits++; 128 129 // if exponent is odd, scale coefficient by 10 130 scale = 67 - digits; 131 exponent_q = exponent_x - scale; 132 scale += (exponent_q & 1); // exp. bias is even 133 134 if (scale > 38) { 135 T128 = power10_table_128[scale - 37]; 136 __mul_128x128_low (CX1, CX, T128); 137 138 TP128 = power10_table_128[37]; 139 __mul_128x128_to_256 (C256, CX1, TP128); 140 } else { 141 T128 = power10_table_128[scale]; 142 __mul_128x128_to_256 (C256, CX, T128); 143 } 144 145 146 // 4*C256 147 C4.w[3] = (C256.w[3] << 2) | (C256.w[2] >> 62); 148 C4.w[2] = (C256.w[2] << 2) | (C256.w[1] >> 62); 149 C4.w[1] = (C256.w[1] << 2) | (C256.w[0] >> 62); 150 C4.w[0] = C256.w[0] << 2; 151 152 long_sqrt128 (&CS, C256); 153 154 #ifndef IEEE_ROUND_NEAREST 155 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY 156 if (!((rnd_mode) & 3)) { 157 #endif 158 #endif 159 // compare to midpoints 160 CSM.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63); 161 CSM.w[0] = (CS.w[0] + CS.w[0]) | 1; 162 // CSM^2 163 //__mul_128x128_to_256(M256, CSM, CSM); 164 __sqr128_to_256 (M256, CSM); 165 166 if (C4.w[3] > M256.w[3] 167 || (C4.w[3] == M256.w[3] 168 && (C4.w[2] > M256.w[2] 169 || (C4.w[2] == M256.w[2] 170 && (C4.w[1] > M256.w[1] 171 || (C4.w[1] == M256.w[1] 172 && C4.w[0] > M256.w[0])))))) { 173 // round up 174 CS.w[0]++; 175 if (!CS.w[0]) 176 CS.w[1]++; 177 } else { 178 C8.w[1] = (CS.w[1] << 3) | (CS.w[0] >> 61); 179 C8.w[0] = CS.w[0] << 3; 180 // M256 - 8*CSM 181 __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]); 182 __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry); 183 __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry); 184 M256.w[3] = M256.w[3] - Carry; 185 186 // if CSM' > C256, round up 187 if (M256.w[3] > C4.w[3] 188 || (M256.w[3] == C4.w[3] 189 && (M256.w[2] > C4.w[2] 190 || (M256.w[2] == C4.w[2] 191 && (M256.w[1] > C4.w[1] 192 || (M256.w[1] == C4.w[1] 193 && M256.w[0] > C4.w[0])))))) { 194 // round down 195 if (!CS.w[0]) 196 CS.w[1]--; 197 CS.w[0]--; 198 } 199 } 200 #ifndef IEEE_ROUND_NEAREST 201 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY 202 } else { 203 __sqr128_to_256 (M256, CS); 204 C8.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63); 205 C8.w[0] = CS.w[0] << 1; 206 if (M256.w[3] > C256.w[3] 207 || (M256.w[3] == C256.w[3] 208 && (M256.w[2] > C256.w[2] 209 || (M256.w[2] == C256.w[2] 210 && (M256.w[1] > C256.w[1] 211 || (M256.w[1] == C256.w[1] 212 && M256.w[0] > C256.w[0])))))) { 213 __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]); 214 __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry); 215 __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry); 216 M256.w[3] = M256.w[3] - Carry; 217 M256.w[0]++; 218 if (!M256.w[0]) { 219 M256.w[1]++; 220 if (!M256.w[1]) { 221 M256.w[2]++; 222 if (!M256.w[2]) 223 M256.w[3]++; 224 } 225 } 226 227 if (!CS.w[0]) 228 CS.w[1]--; 229 CS.w[0]--; 230 231 if (M256.w[3] > C256.w[3] 232 || (M256.w[3] == C256.w[3] 233 && (M256.w[2] > C256.w[2] 234 || (M256.w[2] == C256.w[2] 235 && (M256.w[1] > C256.w[1] 236 || (M256.w[1] == C256.w[1] 237 && M256.w[0] > C256.w[0])))))) { 238 239 if (!CS.w[0]) 240 CS.w[1]--; 241 CS.w[0]--; 242 } 243 } 244 245 else { 246 __add_carry_out (M256.w[0], Carry, M256.w[0], C8.w[0]); 247 __add_carry_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry); 248 __add_carry_in_out (M256.w[2], Carry, M256.w[2], 0, Carry); 249 M256.w[3] = M256.w[3] + Carry; 250 M256.w[0]++; 251 if (!M256.w[0]) { 252 M256.w[1]++; 253 if (!M256.w[1]) { 254 M256.w[2]++; 255 if (!M256.w[2]) 256 M256.w[3]++; 257 } 258 } 259 if (M256.w[3] < C256.w[3] 260 || (M256.w[3] == C256.w[3] 261 && (M256.w[2] < C256.w[2] 262 || (M256.w[2] == C256.w[2] 263 && (M256.w[1] < C256.w[1] 264 || (M256.w[1] == C256.w[1] 265 && M256.w[0] <= C256.w[0])))))) { 266 267 CS.w[0]++; 268 if (!CS.w[0]) 269 CS.w[1]++; 270 } 271 } 272 // RU? 273 if ((rnd_mode) == ROUNDING_UP) { 274 CS.w[0]++; 275 if (!CS.w[0]) 276 CS.w[1]++; 277 } 278 279 } 280 #endif 281 #endif 282 283 #ifdef SET_STATUS_FLAGS 284 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 285 #endif 286 get_BID128_fast (&res, 0, 287 (exponent_q + DECIMAL_EXPONENT_BIAS_128) >> 1, CS); 288 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 289 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 290 #endif 291 BID_RETURN (res); 292 } 293 294 295 296 BID128_FUNCTION_ARGTYPE1 (bid128d_sqrt, UINT64, x) 297 298 UINT256 M256, C256, C4, C8; 299 UINT128 CX, CX1, CX2, A10, S2, T128, TP128, CS, CSM, res; 300 UINT64 sign_x, Carry; 301 SINT64 D; 302 int_float fx, f64; 303 int exponent_x, bin_expon_cx; 304 int digits, scale, exponent_q; 305 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 306 fexcept_t binaryflags = 0; 307 #endif 308 309 // unpack arguments, check for NaN or Infinity 310 // unpack arguments, check for NaN or Infinity 311 CX.w[1] = 0; 312 if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], x)) { 313 res.w[1] = CX.w[0]; 314 res.w[0] = 0; 315 // NaN ? 316 if ((x & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 317 #ifdef SET_STATUS_FLAGS 318 if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN 319 __set_status_flags (pfpsf, INVALID_EXCEPTION); 320 #endif 321 res.w[0] = (CX.w[0] & 0x0003ffffffffffffull); 322 __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]); 323 res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull); 324 BID_RETURN (res); 325 } 326 // x is Infinity? 327 if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) { 328 if (sign_x) { 329 // -Inf, return NaN 330 res.w[1] = 0x7c00000000000000ull; 331 #ifdef SET_STATUS_FLAGS 332 __set_status_flags (pfpsf, INVALID_EXCEPTION); 333 #endif 334 } 335 BID_RETURN (res); 336 } 337 // x is 0 otherwise 338 339 exponent_x = 340 exponent_x - DECIMAL_EXPONENT_BIAS + DECIMAL_EXPONENT_BIAS_128; 341 res.w[1] = 342 sign_x | ((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1) 343 << 49); 344 res.w[0] = 0; 345 BID_RETURN (res); 346 } 347 if (sign_x) { 348 res.w[1] = 0x7c00000000000000ull; 349 res.w[0] = 0; 350 #ifdef SET_STATUS_FLAGS 351 __set_status_flags (pfpsf, INVALID_EXCEPTION); 352 #endif 353 BID_RETURN (res); 354 } 355 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 356 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS); 357 #endif 358 exponent_x = 359 exponent_x - DECIMAL_EXPONENT_BIAS + DECIMAL_EXPONENT_BIAS_128; 360 361 // 2^64 362 f64.i = 0x5f800000; 363 364 // fx ~ CX 365 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0]; 366 bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f; 367 digits = estimate_decimal_digits[bin_expon_cx]; 368 369 A10 = CX; 370 if (exponent_x & 1) { 371 A10.w[1] = (CX.w[1] << 3) | (CX.w[0] >> 61); 372 A10.w[0] = CX.w[0] << 3; 373 CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63); 374 CX2.w[0] = CX.w[0] << 1; 375 __add_128_128 (A10, A10, CX2); 376 } 377 378 CS.w[0] = short_sqrt128 (A10); 379 CS.w[1] = 0; 380 // check for exact result 381 if (CS.w[0] * CS.w[0] == A10.w[0]) { 382 __mul_64x64_to_128_fast (S2, CS.w[0], CS.w[0]); 383 if (S2.w[1] == A10.w[1]) { 384 get_BID128_very_fast (&res, 0, 385 (exponent_x + DECIMAL_EXPONENT_BIAS_128) >> 1, 386 CS); 387 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 388 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 389 #endif 390 BID_RETURN (res); 391 } 392 } 393 // get number of digits in CX 394 D = CX.w[1] - power10_index_binexp_128[bin_expon_cx].w[1]; 395 if (D > 0 396 || (!D && CX.w[0] >= power10_index_binexp_128[bin_expon_cx].w[0])) 397 digits++; 398 399 // if exponent is odd, scale coefficient by 10 400 scale = 67 - digits; 401 exponent_q = exponent_x - scale; 402 scale += (exponent_q & 1); // exp. bias is even 403 404 if (scale > 38) { 405 T128 = power10_table_128[scale - 37]; 406 __mul_128x128_low (CX1, CX, T128); 407 408 TP128 = power10_table_128[37]; 409 __mul_128x128_to_256 (C256, CX1, TP128); 410 } else { 411 T128 = power10_table_128[scale]; 412 __mul_128x128_to_256 (C256, CX, T128); 413 } 414 415 416 // 4*C256 417 C4.w[3] = (C256.w[3] << 2) | (C256.w[2] >> 62); 418 C4.w[2] = (C256.w[2] << 2) | (C256.w[1] >> 62); 419 C4.w[1] = (C256.w[1] << 2) | (C256.w[0] >> 62); 420 C4.w[0] = C256.w[0] << 2; 421 422 long_sqrt128 (&CS, C256); 423 424 #ifndef IEEE_ROUND_NEAREST 425 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY 426 if (!((rnd_mode) & 3)) { 427 #endif 428 #endif 429 // compare to midpoints 430 CSM.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63); 431 CSM.w[0] = (CS.w[0] + CS.w[0]) | 1; 432 // CSM^2 433 //__mul_128x128_to_256(M256, CSM, CSM); 434 __sqr128_to_256 (M256, CSM); 435 436 if (C4.w[3] > M256.w[3] 437 || (C4.w[3] == M256.w[3] 438 && (C4.w[2] > M256.w[2] 439 || (C4.w[2] == M256.w[2] 440 && (C4.w[1] > M256.w[1] 441 || (C4.w[1] == M256.w[1] 442 && C4.w[0] > M256.w[0])))))) { 443 // round up 444 CS.w[0]++; 445 if (!CS.w[0]) 446 CS.w[1]++; 447 } else { 448 C8.w[1] = (CS.w[1] << 3) | (CS.w[0] >> 61); 449 C8.w[0] = CS.w[0] << 3; 450 // M256 - 8*CSM 451 __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]); 452 __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry); 453 __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry); 454 M256.w[3] = M256.w[3] - Carry; 455 456 // if CSM' > C256, round up 457 if (M256.w[3] > C4.w[3] 458 || (M256.w[3] == C4.w[3] 459 && (M256.w[2] > C4.w[2] 460 || (M256.w[2] == C4.w[2] 461 && (M256.w[1] > C4.w[1] 462 || (M256.w[1] == C4.w[1] 463 && M256.w[0] > C4.w[0])))))) { 464 // round down 465 if (!CS.w[0]) 466 CS.w[1]--; 467 CS.w[0]--; 468 } 469 } 470 #ifndef IEEE_ROUND_NEAREST 471 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY 472 } else { 473 __sqr128_to_256 (M256, CS); 474 C8.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63); 475 C8.w[0] = CS.w[0] << 1; 476 if (M256.w[3] > C256.w[3] 477 || (M256.w[3] == C256.w[3] 478 && (M256.w[2] > C256.w[2] 479 || (M256.w[2] == C256.w[2] 480 && (M256.w[1] > C256.w[1] 481 || (M256.w[1] == C256.w[1] 482 && M256.w[0] > C256.w[0])))))) { 483 __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]); 484 __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry); 485 __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry); 486 M256.w[3] = M256.w[3] - Carry; 487 M256.w[0]++; 488 if (!M256.w[0]) { 489 M256.w[1]++; 490 if (!M256.w[1]) { 491 M256.w[2]++; 492 if (!M256.w[2]) 493 M256.w[3]++; 494 } 495 } 496 497 if (!CS.w[0]) 498 CS.w[1]--; 499 CS.w[0]--; 500 501 if (M256.w[3] > C256.w[3] 502 || (M256.w[3] == C256.w[3] 503 && (M256.w[2] > C256.w[2] 504 || (M256.w[2] == C256.w[2] 505 && (M256.w[1] > C256.w[1] 506 || (M256.w[1] == C256.w[1] 507 && M256.w[0] > C256.w[0])))))) { 508 509 if (!CS.w[0]) 510 CS.w[1]--; 511 CS.w[0]--; 512 } 513 } 514 515 else { 516 __add_carry_out (M256.w[0], Carry, M256.w[0], C8.w[0]); 517 __add_carry_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry); 518 __add_carry_in_out (M256.w[2], Carry, M256.w[2], 0, Carry); 519 M256.w[3] = M256.w[3] + Carry; 520 M256.w[0]++; 521 if (!M256.w[0]) { 522 M256.w[1]++; 523 if (!M256.w[1]) { 524 M256.w[2]++; 525 if (!M256.w[2]) 526 M256.w[3]++; 527 } 528 } 529 if (M256.w[3] < C256.w[3] 530 || (M256.w[3] == C256.w[3] 531 && (M256.w[2] < C256.w[2] 532 || (M256.w[2] == C256.w[2] 533 && (M256.w[1] < C256.w[1] 534 || (M256.w[1] == C256.w[1] 535 && M256.w[0] <= C256.w[0])))))) { 536 537 CS.w[0]++; 538 if (!CS.w[0]) 539 CS.w[1]++; 540 } 541 } 542 // RU? 543 if ((rnd_mode) == ROUNDING_UP) { 544 CS.w[0]++; 545 if (!CS.w[0]) 546 CS.w[1]++; 547 } 548 549 } 550 #endif 551 #endif 552 553 #ifdef SET_STATUS_FLAGS 554 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 555 #endif 556 get_BID128_fast (&res, 0, (exponent_q + DECIMAL_EXPONENT_BIAS_128) >> 1, 557 CS); 558 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 559 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 560 #endif 561 BID_RETURN (res); 562 563 564 } 565