1 //===-- include/flang/Evaluate/integer.h ------------------------*- C++ -*-===// 2 // 3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. 4 // See https://llvm.org/LICENSE.txt for license information. 5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception 6 // 7 //===----------------------------------------------------------------------===// 8 9 #ifndef FORTRAN_EVALUATE_INTEGER_H_ 10 #define FORTRAN_EVALUATE_INTEGER_H_ 11 12 // Emulates binary integers of an arbitrary (but fixed) bit size for use 13 // when the host C++ environment does not support that size or when the 14 // full suite of Fortran's integer intrinsic scalar functions are needed. 15 // The data model is typeless, so signed* and unsigned operations 16 // are distinguished from each other with distinct member function interfaces. 17 // (*"Signed" here means two's-complement, just to be clear. Ones'-complement 18 // and signed-magnitude encodings appear to be extinct in 2018.) 19 20 #include "flang/Common/bit-population-count.h" 21 #include "flang/Common/leading-zero-bit-count.h" 22 #include "flang/Evaluate/common.h" 23 #include <cinttypes> 24 #include <climits> 25 #include <cstddef> 26 #include <cstdint> 27 #include <string> 28 #include <type_traits> 29 30 // Some environments, viz. clang on Darwin, allow the macro HUGE 31 // to leak out of <math.h> even when it is never directly included. 32 #undef HUGE 33 34 namespace Fortran::evaluate::value { 35 36 // Implements an integer as an assembly of smaller host integer parts 37 // that constitute the digits of a large-radix fixed-point number. 38 // For best performance, the type of these parts should be half of the 39 // size of the largest efficient integer supported by the host processor. 40 // These parts are stored in either little- or big-endian order, which can 41 // match that of the host's endianness or not; but if the ordering matches 42 // that of the host, raw host data can be overlaid with a properly configured 43 // instance of this class and used in situ. 44 // To facilitate exhaustive testing of what would otherwise be more rare 45 // edge cases, this class template may be configured to use other part 46 // types &/or partial fields in the parts. The radix (i.e., the number 47 // of possible values in a part), however, must be a power of two; this 48 // template class is not generalized to enable, say, decimal arithmetic. 49 // Member functions that correspond to Fortran intrinsic functions are 50 // named accordingly in ALL CAPS so that they can be referenced easily in 51 // the language standard. 52 template <int BITS, bool IS_LITTLE_ENDIAN = isHostLittleEndian, 53 int PARTBITS = BITS <= 32 ? BITS : 32, 54 typename PART = HostUnsignedInt<PARTBITS>, 55 typename BIGPART = HostUnsignedInt<PARTBITS * 2>> 56 class Integer { 57 public: 58 static constexpr int bits{BITS}; 59 static constexpr int partBits{PARTBITS}; 60 using Part = PART; 61 using BigPart = BIGPART; 62 static_assert(std::is_integral_v<Part>); 63 static_assert(std::is_unsigned_v<Part>); 64 static_assert(std::is_integral_v<BigPart>); 65 static_assert(std::is_unsigned_v<BigPart>); 66 static_assert(CHAR_BIT * sizeof(BigPart) >= 2 * partBits); 67 static constexpr bool littleEndian{IS_LITTLE_ENDIAN}; 68 69 private: 70 static constexpr int maxPartBits{CHAR_BIT * sizeof(Part)}; 71 static_assert(partBits > 0 && partBits <= maxPartBits); 72 static constexpr int extraPartBits{maxPartBits - partBits}; 73 static constexpr int parts{(bits + partBits - 1) / partBits}; 74 static_assert(parts >= 1); 75 static constexpr int extraTopPartBits{ 76 extraPartBits + (parts * partBits) - bits}; 77 static constexpr int topPartBits{maxPartBits - extraTopPartBits}; 78 static_assert(topPartBits > 0 && topPartBits <= partBits); 79 static_assert((parts - 1) * partBits + topPartBits == bits); 80 static constexpr Part partMask{static_cast<Part>(~0) >> extraPartBits}; 81 static constexpr Part topPartMask{static_cast<Part>(~0) >> extraTopPartBits}; 82 83 public: 84 // Some types used for member function results 85 struct ValueWithOverflow { 86 Integer value; 87 bool overflow; 88 }; 89 90 struct ValueWithCarry { 91 Integer value; 92 bool carry; 93 }; 94 95 struct Product { SignedMultiplicationOverflowedProduct96 bool SignedMultiplicationOverflowed() const { 97 return lower.IsNegative() ? (upper.POPCNT() != bits) : !upper.IsZero(); 98 } 99 Integer upper, lower; 100 }; 101 102 struct QuotientWithRemainder { 103 Integer quotient, remainder; 104 bool divisionByZero, overflow; 105 }; 106 107 struct PowerWithErrors { 108 Integer power; 109 bool divisionByZero{false}, overflow{false}, zeroToZero{false}; 110 }; 111 112 // Constructors and value-generating static functions Integer()113 constexpr Integer() { Clear(); } // default constructor: zero 114 constexpr Integer(const Integer &) = default; 115 constexpr Integer(Integer &&) = default; 116 117 // C++'s integral types can all be converted to Integer 118 // with silent truncation. 119 template <typename INT, typename = std::enable_if_t<std::is_integral_v<INT>>> Integer(INT n)120 constexpr Integer(INT n) { 121 constexpr int nBits = CHAR_BIT * sizeof n; 122 if constexpr (nBits < partBits) { 123 if constexpr (std::is_unsigned_v<INT>) { 124 // Zero-extend an unsigned smaller value. 125 SetLEPart(0, n); 126 for (int j{1}; j < parts; ++j) { 127 SetLEPart(j, 0); 128 } 129 } else { 130 // n has a signed type smaller than the usable 131 // bits in a Part. 132 // Avoid conversions that change both size and sign. 133 using SignedPart = std::make_signed_t<Part>; 134 Part p = static_cast<SignedPart>(n); 135 SetLEPart(0, p); 136 if constexpr (parts > 1) { 137 Part signExtension = static_cast<SignedPart>(-(n < 0)); 138 for (int j{1}; j < parts; ++j) { 139 SetLEPart(j, signExtension); 140 } 141 } 142 } 143 } else { 144 // n has some integral type no smaller than the usable 145 // bits in a Part. 146 // Ensure that all shifts are smaller than a whole word. 147 if constexpr (std::is_unsigned_v<INT>) { 148 for (int j{0}; j < parts; ++j) { 149 SetLEPart(j, static_cast<Part>(n)); 150 if constexpr (nBits > partBits) { 151 n >>= partBits; 152 } else { 153 n = 0; 154 } 155 } 156 } else { 157 INT signExtension{-(n < 0)}; 158 static_assert(nBits >= partBits); 159 if constexpr (nBits > partBits) { 160 signExtension <<= nBits - partBits; 161 for (int j{0}; j < parts; ++j) { 162 SetLEPart(j, static_cast<Part>(n)); 163 n >>= partBits; 164 n |= signExtension; 165 } 166 } else { 167 SetLEPart(0, static_cast<Part>(n)); 168 for (int j{1}; j < parts; ++j) { 169 SetLEPart(j, static_cast<Part>(signExtension)); 170 } 171 } 172 } 173 } 174 } 175 176 constexpr Integer &operator=(const Integer &) = default; 177 178 constexpr bool operator<(const Integer &that) const { 179 return CompareSigned(that) == Ordering::Less; 180 } 181 constexpr bool operator<=(const Integer &that) const { 182 return CompareSigned(that) != Ordering::Greater; 183 } 184 constexpr bool operator==(const Integer &that) const { 185 return CompareSigned(that) == Ordering::Equal; 186 } 187 constexpr bool operator!=(const Integer &that) const { 188 return !(*this == that); 189 } 190 constexpr bool operator>=(const Integer &that) const { 191 return CompareSigned(that) != Ordering::Less; 192 } 193 constexpr bool operator>(const Integer &that) const { 194 return CompareSigned(that) == Ordering::Greater; 195 } 196 197 // Left-justified mask (e.g., MASKL(1) has only its sign bit set) MASKL(int places)198 static constexpr Integer MASKL(int places) { 199 if (places <= 0) { 200 return {}; 201 } else if (places >= bits) { 202 return MASKR(bits); 203 } else { 204 return MASKR(bits - places).NOT(); 205 } 206 } 207 208 // Right-justified mask (e.g., MASKR(1) == 1, MASKR(2) == 3, &c.) MASKR(int places)209 static constexpr Integer MASKR(int places) { 210 Integer result{nullptr}; 211 int j{0}; 212 for (; j + 1 < parts && places >= partBits; ++j, places -= partBits) { 213 result.LEPart(j) = partMask; 214 } 215 if (places > 0) { 216 if (j + 1 < parts) { 217 result.LEPart(j++) = partMask >> (partBits - places); 218 } else if (j + 1 == parts) { 219 if (places >= topPartBits) { 220 result.LEPart(j++) = topPartMask; 221 } else { 222 result.LEPart(j++) = topPartMask >> (topPartBits - places); 223 } 224 } 225 } 226 for (; j < parts; ++j) { 227 result.LEPart(j) = 0; 228 } 229 return result; 230 } 231 232 static constexpr ValueWithOverflow Read( 233 const char *&pp, std::uint64_t base = 10, bool isSigned = false) { 234 Integer result; 235 bool overflow{false}; 236 const char *p{pp}; 237 while (*p == ' ' || *p == '\t') { 238 ++p; 239 } 240 bool negate{*p == '-'}; 241 if (negate || *p == '+') { 242 while (*++p == ' ' || *p == '\t') { 243 } 244 } 245 Integer radix{base}; 246 // This code makes assumptions about local contiguity in regions of the 247 // character set and only works up to base 36. These assumptions hold 248 // for all current combinations of surviving character sets (ASCII, UTF-8, 249 // EBCDIC) and the bases used in Fortran source and formatted I/O 250 // (viz., 2, 8, 10, & 16). But: management thought that a disclaimer 251 // might be needed here to warn future users of this code about these 252 // assumptions, so here you go, future programmer in some postapocalyptic 253 // hellscape, and best of luck with the inexorable killer robots. 254 for (; std::uint64_t digit = *p; ++p) { 255 if (digit >= '0' && digit <= '9' && digit < '0' + base) { 256 digit -= '0'; 257 } else if (base > 10 && digit >= 'A' && digit < 'A' + base - 10) { 258 digit -= 'A' - 10; 259 } else if (base > 10 && digit >= 'a' && digit < 'a' + base - 10) { 260 digit -= 'a' - 10; 261 } else { 262 break; 263 } 264 Product shifted{result.MultiplyUnsigned(radix)}; 265 overflow |= !shifted.upper.IsZero(); 266 ValueWithCarry next{shifted.lower.AddUnsigned(Integer{digit})}; 267 overflow |= next.carry; 268 result = next.value; 269 } 270 pp = p; 271 if (negate) { 272 result = result.Negate().value; 273 overflow |= isSigned && !result.IsNegative() && !result.IsZero(); 274 } else { 275 overflow |= isSigned && result.IsNegative(); 276 } 277 return {result, overflow}; 278 } 279 280 template <typename FROM> ConvertUnsigned(const FROM & that)281 static constexpr ValueWithOverflow ConvertUnsigned(const FROM &that) { 282 std::uint64_t field{that.ToUInt64()}; 283 ValueWithOverflow result{field, false}; 284 if constexpr (bits < 64) { 285 result.overflow = (field >> bits) != 0; 286 } 287 for (int j{64}; j < that.bits && !result.overflow; j += 64) { 288 field = that.SHIFTR(j).ToUInt64(); 289 if (bits <= j) { 290 result.overflow = field != 0; 291 } else { 292 result.value = result.value.IOR(Integer{field}.SHIFTL(j)); 293 if (bits < j + 64) { 294 result.overflow = (field >> (bits - j)) != 0; 295 } 296 } 297 } 298 return result; 299 } 300 301 template <typename FROM> ConvertSigned(const FROM & that)302 static constexpr ValueWithOverflow ConvertSigned(const FROM &that) { 303 ValueWithOverflow result{ConvertUnsigned(that)}; 304 if constexpr (bits > FROM::bits) { 305 if (that.IsNegative()) { 306 result.value = result.value.IOR(MASKL(bits - FROM::bits)); 307 } 308 result.overflow = false; 309 } else if constexpr (bits < FROM::bits) { 310 auto back{FROM::template ConvertSigned(result.value)}; 311 result.overflow = back.value.CompareUnsigned(that) != Ordering::Equal; 312 } 313 return result; 314 } 315 UnsignedDecimal()316 std::string UnsignedDecimal() const { 317 if constexpr (bits < 4) { 318 char digit = '0' + ToUInt64(); 319 return {digit}; 320 } else if (IsZero()) { 321 return {'0'}; 322 } else { 323 QuotientWithRemainder qr{DivideUnsigned(10)}; 324 char digit = '0' + qr.remainder.ToUInt64(); 325 if (qr.quotient.IsZero()) { 326 return {digit}; 327 } else { 328 return qr.quotient.UnsignedDecimal() + digit; 329 } 330 } 331 } 332 SignedDecimal()333 std::string SignedDecimal() const { 334 if (IsNegative()) { 335 return std::string{'-'} + Negate().value.UnsignedDecimal(); 336 } else { 337 return UnsignedDecimal(); 338 } 339 } 340 341 // Omits a leading "0x". Hexadecimal()342 std::string Hexadecimal() const { 343 std::string result; 344 int digits{(bits + 3) >> 2}; 345 for (int j{0}; j < digits; ++j) { 346 int pos{(digits - 1 - j) * 4}; 347 char nybble = IBITS(pos, 4).ToUInt64(); 348 if (nybble != 0 || !result.empty() || j + 1 == digits) { 349 char digit = '0' + nybble; 350 if (digit > '9') { 351 digit += 'a' - ('9' + 1); 352 } 353 result += digit; 354 } 355 } 356 return result; 357 } 358 359 static constexpr int DIGITS{bits - 1}; // don't count the sign bit HUGE()360 static constexpr Integer HUGE() { return MASKR(bits - 1); } Least()361 static constexpr Integer Least() { return MASKL(1); } 362 static constexpr int RANGE{// in the sense of SELECTED_INT_KIND 363 // This magic value is LOG10(2.)*1E12. 364 static_cast<int>(((bits - 1) * 301029995664) / 1000000000000)}; 365 IsZero()366 constexpr bool IsZero() const { 367 for (int j{0}; j < parts; ++j) { 368 if (part_[j] != 0) { 369 return false; 370 } 371 } 372 return true; 373 } 374 IsNegative()375 constexpr bool IsNegative() const { 376 return (LEPart(parts - 1) >> (topPartBits - 1)) & 1; 377 } 378 CompareToZeroSigned()379 constexpr Ordering CompareToZeroSigned() const { 380 if (IsNegative()) { 381 return Ordering::Less; 382 } else if (IsZero()) { 383 return Ordering::Equal; 384 } else { 385 return Ordering::Greater; 386 } 387 } 388 389 // Count the number of contiguous most-significant bit positions 390 // that are clear. LEADZ()391 constexpr int LEADZ() const { 392 if (LEPart(parts - 1) != 0) { 393 int lzbc{common::LeadingZeroBitCount(LEPart(parts - 1))}; 394 return lzbc - extraTopPartBits; 395 } 396 int upperZeroes{topPartBits}; 397 for (int j{1}; j < parts; ++j) { 398 if (Part p{LEPart(parts - 1 - j)}) { 399 int lzbc{common::LeadingZeroBitCount(p)}; 400 return upperZeroes + lzbc - extraPartBits; 401 } 402 upperZeroes += partBits; 403 } 404 return bits; 405 } 406 407 // Count the number of bit positions that are set. POPCNT()408 constexpr int POPCNT() const { 409 int count{0}; 410 for (int j{0}; j < parts; ++j) { 411 count += common::BitPopulationCount(part_[j]); 412 } 413 return count; 414 } 415 416 // True when POPCNT is odd. POPPAR()417 constexpr bool POPPAR() const { return POPCNT() & 1; } 418 TRAILZ()419 constexpr int TRAILZ() const { 420 auto minus1{AddUnsigned(MASKR(bits))}; // { x-1, carry = x > 0 } 421 if (!minus1.carry) { 422 return bits; // was zero 423 } else { 424 // x ^ (x-1) has all bits set at and below original least-order set bit. 425 return IEOR(minus1.value).POPCNT() - 1; 426 } 427 } 428 BTEST(int pos)429 constexpr bool BTEST(int pos) const { 430 if (pos < 0 || pos >= bits) { 431 return false; 432 } else { 433 return (LEPart(pos / partBits) >> (pos % partBits)) & 1; 434 } 435 } 436 CompareUnsigned(const Integer & y)437 constexpr Ordering CompareUnsigned(const Integer &y) const { 438 for (int j{parts}; j-- > 0;) { 439 if (LEPart(j) > y.LEPart(j)) { 440 return Ordering::Greater; 441 } 442 if (LEPart(j) < y.LEPart(j)) { 443 return Ordering::Less; 444 } 445 } 446 return Ordering::Equal; 447 } 448 BGE(const Integer & y)449 constexpr bool BGE(const Integer &y) const { 450 return CompareUnsigned(y) != Ordering::Less; 451 } BGT(const Integer & y)452 constexpr bool BGT(const Integer &y) const { 453 return CompareUnsigned(y) == Ordering::Greater; 454 } BLE(const Integer & y)455 constexpr bool BLE(const Integer &y) const { return !BGT(y); } BLT(const Integer & y)456 constexpr bool BLT(const Integer &y) const { return !BGE(y); } 457 CompareSigned(const Integer & y)458 constexpr Ordering CompareSigned(const Integer &y) const { 459 bool isNegative{IsNegative()}; 460 if (isNegative != y.IsNegative()) { 461 return isNegative ? Ordering::Less : Ordering::Greater; 462 } 463 return CompareUnsigned(y); 464 } 465 ToUInt()466 template <typename UINT = std::uint64_t> constexpr UINT ToUInt() const { 467 UINT n{LEPart(0)}; 468 std::size_t filled{partBits}; 469 constexpr std::size_t maxBits{CHAR_BIT * sizeof n}; 470 for (int j{1}; filled < maxBits && j < parts; ++j, filled += partBits) { 471 n |= UINT{LEPart(j)} << filled; 472 } 473 return n; 474 } 475 476 template <typename SINT = std::int64_t, typename UINT = std::uint64_t> ToSInt()477 constexpr SINT ToSInt() const { 478 SINT n = ToUInt<UINT>(); 479 constexpr std::size_t maxBits{CHAR_BIT * sizeof n}; 480 if constexpr (bits < maxBits) { 481 n |= -(n >> (bits - 1)) << bits; 482 } 483 return n; 484 } 485 ToUInt64()486 constexpr std::uint64_t ToUInt64() const { return ToUInt<std::uint64_t>(); } 487 ToInt64()488 constexpr std::int64_t ToInt64() const { 489 return ToSInt<std::int64_t, std::uint64_t>(); 490 } 491 492 // Ones'-complement (i.e., C's ~) NOT()493 constexpr Integer NOT() const { 494 Integer result{nullptr}; 495 for (int j{0}; j < parts; ++j) { 496 result.SetLEPart(j, ~LEPart(j)); 497 } 498 return result; 499 } 500 501 // Two's-complement negation (-x = ~x + 1). 502 // An overflow flag accompanies the result, and will be true when the 503 // operand is the most negative signed number (MASKL(1)). Negate()504 constexpr ValueWithOverflow Negate() const { 505 Integer result{nullptr}; 506 Part carry{1}; 507 for (int j{0}; j + 1 < parts; ++j) { 508 Part newCarry{LEPart(j) == 0 && carry}; 509 result.SetLEPart(j, ~LEPart(j) + carry); 510 carry = newCarry; 511 } 512 Part top{LEPart(parts - 1)}; 513 result.SetLEPart(parts - 1, ~top + carry); 514 bool overflow{top != 0 && result.LEPart(parts - 1) == top}; 515 return {result, overflow}; 516 } 517 ABS()518 constexpr ValueWithOverflow ABS() const { 519 if (IsNegative()) { 520 return Negate(); 521 } else { 522 return {*this, false}; 523 } 524 } 525 526 // Shifts the operand left when the count is positive, right when negative. 527 // Vacated bit positions are filled with zeroes. ISHFT(int count)528 constexpr Integer ISHFT(int count) const { 529 if (count < 0) { 530 return SHIFTR(-count); 531 } else { 532 return SHIFTL(count); 533 } 534 } 535 536 // Left shift with zero fill. SHIFTL(int count)537 constexpr Integer SHIFTL(int count) const { 538 if (count <= 0) { 539 return *this; 540 } else { 541 Integer result{nullptr}; 542 int shiftParts{count / partBits}; 543 int bitShift{count - partBits * shiftParts}; 544 int j{parts - 1}; 545 if (bitShift == 0) { 546 for (; j >= shiftParts; --j) { 547 result.SetLEPart(j, LEPart(j - shiftParts)); 548 } 549 for (; j >= 0; --j) { 550 result.LEPart(j) = 0; 551 } 552 } else { 553 for (; j > shiftParts; --j) { 554 result.SetLEPart(j, 555 ((LEPart(j - shiftParts) << bitShift) | 556 (LEPart(j - shiftParts - 1) >> (partBits - bitShift)))); 557 } 558 if (j == shiftParts) { 559 result.SetLEPart(j, LEPart(0) << bitShift); 560 --j; 561 } 562 for (; j >= 0; --j) { 563 result.LEPart(j) = 0; 564 } 565 } 566 return result; 567 } 568 } 569 570 // Circular shift of a field of least-significant bits. The least-order 571 // "size" bits are shifted circularly in place by "count" positions; 572 // the shift is leftward if count is nonnegative, rightward otherwise. 573 // Higher-order bits are unchanged. 574 constexpr Integer ISHFTC(int count, int size = bits) const { 575 if (count == 0 || size <= 0) { 576 return *this; 577 } 578 if (size > bits) { 579 size = bits; 580 } 581 count %= size; 582 if (count == 0) { 583 return *this; 584 } 585 int middleBits{size - count}, leastBits{count}; 586 if (count < 0) { 587 middleBits = -count; 588 leastBits = size + count; 589 } 590 if (size == bits) { 591 return SHIFTL(leastBits).IOR(SHIFTR(middleBits)); 592 } 593 Integer unchanged{IAND(MASKL(bits - size))}; 594 Integer middle{IAND(MASKR(middleBits)).SHIFTL(leastBits)}; 595 Integer least{SHIFTR(middleBits).IAND(MASKR(leastBits))}; 596 return unchanged.IOR(middle).IOR(least); 597 } 598 599 // Double shifts, aka shifts with specific fill. SHIFTLWithFill(const Integer & fill,int count)600 constexpr Integer SHIFTLWithFill(const Integer &fill, int count) const { 601 if (count <= 0) { 602 return *this; 603 } else if (count >= 2 * bits) { 604 return {}; 605 } else if (count > bits) { 606 return fill.SHIFTL(count - bits); 607 } else if (count == bits) { 608 return fill; 609 } else { 610 return SHIFTL(count).IOR(fill.SHIFTR(bits - count)); 611 } 612 } 613 SHIFTRWithFill(const Integer & fill,int count)614 constexpr Integer SHIFTRWithFill(const Integer &fill, int count) const { 615 if (count <= 0) { 616 return *this; 617 } else if (count >= 2 * bits) { 618 return {}; 619 } else if (count > bits) { 620 return fill.SHIFTR(count - bits); 621 } else if (count == bits) { 622 return fill; 623 } else { 624 return SHIFTR(count).IOR(fill.SHIFTL(bits - count)); 625 } 626 } 627 DSHIFTL(const Integer & fill,int count)628 constexpr Integer DSHIFTL(const Integer &fill, int count) const { 629 // DSHIFTL(I,J) shifts I:J left; the second argument is the right fill. 630 return SHIFTLWithFill(fill, count); 631 } 632 DSHIFTR(const Integer & value,int count)633 constexpr Integer DSHIFTR(const Integer &value, int count) const { 634 // DSHIFTR(I,J) shifts I:J right; the *first* argument is the left fill. 635 return value.SHIFTRWithFill(*this, count); 636 } 637 638 // Vacated upper bits are filled with zeroes. SHIFTR(int count)639 constexpr Integer SHIFTR(int count) const { 640 if (count <= 0) { 641 return *this; 642 } else { 643 Integer result{nullptr}; 644 int shiftParts{count / partBits}; 645 int bitShift{count - partBits * shiftParts}; 646 int j{0}; 647 if (bitShift == 0) { 648 for (; j + shiftParts < parts; ++j) { 649 result.LEPart(j) = LEPart(j + shiftParts); 650 } 651 for (; j < parts; ++j) { 652 result.LEPart(j) = 0; 653 } 654 } else { 655 for (; j + shiftParts + 1 < parts; ++j) { 656 result.SetLEPart(j, 657 (LEPart(j + shiftParts) >> bitShift) | 658 (LEPart(j + shiftParts + 1) << (partBits - bitShift))); 659 } 660 if (j + shiftParts + 1 == parts) { 661 result.LEPart(j++) = LEPart(parts - 1) >> bitShift; 662 } 663 for (; j < parts; ++j) { 664 result.LEPart(j) = 0; 665 } 666 } 667 return result; 668 } 669 } 670 671 // Be advised, an arithmetic (sign-filling) right shift is not 672 // the same as a division by a power of two in all cases. SHIFTA(int count)673 constexpr Integer SHIFTA(int count) const { 674 if (count <= 0) { 675 return *this; 676 } else if (IsNegative()) { 677 return SHIFTR(count).IOR(MASKL(count)); 678 } else { 679 return SHIFTR(count); 680 } 681 } 682 683 // Clears a single bit. IBCLR(int pos)684 constexpr Integer IBCLR(int pos) const { 685 if (pos < 0 || pos >= bits) { 686 return *this; 687 } else { 688 Integer result{*this}; 689 result.LEPart(pos / partBits) &= ~(Part{1} << (pos % partBits)); 690 return result; 691 } 692 } 693 694 // Sets a single bit. IBSET(int pos)695 constexpr Integer IBSET(int pos) const { 696 if (pos < 0 || pos >= bits) { 697 return *this; 698 } else { 699 Integer result{*this}; 700 result.LEPart(pos / partBits) |= Part{1} << (pos % partBits); 701 return result; 702 } 703 } 704 705 // Extracts a field. IBITS(int pos,int size)706 constexpr Integer IBITS(int pos, int size) const { 707 return SHIFTR(pos).IAND(MASKR(size)); 708 } 709 IAND(const Integer & y)710 constexpr Integer IAND(const Integer &y) const { 711 Integer result{nullptr}; 712 for (int j{0}; j < parts; ++j) { 713 result.LEPart(j) = LEPart(j) & y.LEPart(j); 714 } 715 return result; 716 } 717 IOR(const Integer & y)718 constexpr Integer IOR(const Integer &y) const { 719 Integer result{nullptr}; 720 for (int j{0}; j < parts; ++j) { 721 result.LEPart(j) = LEPart(j) | y.LEPart(j); 722 } 723 return result; 724 } 725 IEOR(const Integer & y)726 constexpr Integer IEOR(const Integer &y) const { 727 Integer result{nullptr}; 728 for (int j{0}; j < parts; ++j) { 729 result.LEPart(j) = LEPart(j) ^ y.LEPart(j); 730 } 731 return result; 732 } 733 MERGE_BITS(const Integer & y,const Integer & mask)734 constexpr Integer MERGE_BITS(const Integer &y, const Integer &mask) const { 735 return IAND(mask).IOR(y.IAND(mask.NOT())); 736 } 737 MAX(const Integer & y)738 constexpr Integer MAX(const Integer &y) const { 739 if (CompareSigned(y) == Ordering::Less) { 740 return y; 741 } else { 742 return *this; 743 } 744 } 745 MIN(const Integer & y)746 constexpr Integer MIN(const Integer &y) const { 747 if (CompareSigned(y) == Ordering::Less) { 748 return *this; 749 } else { 750 return y; 751 } 752 } 753 754 // Unsigned addition with carry. 755 constexpr ValueWithCarry AddUnsigned( 756 const Integer &y, bool carryIn = false) const { 757 Integer sum{nullptr}; 758 BigPart carry{carryIn}; 759 for (int j{0}; j + 1 < parts; ++j) { 760 carry += LEPart(j); 761 carry += y.LEPart(j); 762 sum.SetLEPart(j, carry); 763 carry >>= partBits; 764 } 765 carry += LEPart(parts - 1); 766 carry += y.LEPart(parts - 1); 767 sum.SetLEPart(parts - 1, carry); 768 return {sum, carry > topPartMask}; 769 } 770 AddSigned(const Integer & y)771 constexpr ValueWithOverflow AddSigned(const Integer &y) const { 772 bool isNegative{IsNegative()}; 773 bool sameSign{isNegative == y.IsNegative()}; 774 ValueWithCarry sum{AddUnsigned(y)}; 775 bool overflow{sameSign && sum.value.IsNegative() != isNegative}; 776 return {sum.value, overflow}; 777 } 778 SubtractSigned(const Integer & y)779 constexpr ValueWithOverflow SubtractSigned(const Integer &y) const { 780 bool isNegative{IsNegative()}; 781 bool sameSign{isNegative == y.IsNegative()}; 782 ValueWithCarry diff{AddUnsigned(y.Negate().value)}; 783 bool overflow{!sameSign && diff.value.IsNegative() != isNegative}; 784 return {diff.value, overflow}; 785 } 786 787 // MAX(X-Y, 0) DIM(const Integer & y)788 constexpr Integer DIM(const Integer &y) const { 789 if (CompareSigned(y) != Ordering::Greater) { 790 return {}; 791 } else { 792 return SubtractSigned(y).value; 793 } 794 } 795 SIGN(bool toNegative)796 constexpr ValueWithOverflow SIGN(bool toNegative) const { 797 if (toNegative == IsNegative()) { 798 return {*this, false}; 799 } else if (toNegative) { 800 return Negate(); 801 } else { 802 return ABS(); 803 } 804 } 805 SIGN(const Integer & sign)806 constexpr ValueWithOverflow SIGN(const Integer &sign) const { 807 return SIGN(sign.IsNegative()); 808 } 809 MultiplyUnsigned(const Integer & y)810 constexpr Product MultiplyUnsigned(const Integer &y) const { 811 Part product[2 * parts]{}; // little-endian full product 812 for (int j{0}; j < parts; ++j) { 813 if (Part xpart{LEPart(j)}) { 814 for (int k{0}; k < parts; ++k) { 815 if (Part ypart{y.LEPart(k)}) { 816 BigPart xy{xpart}; 817 xy *= ypart; 818 #if defined __GNUC__ && __GNUC__ < 8 819 // && to < (2 * parts) was added to avoid GCC < 8 build failure on 820 // -Werror=array-bounds. This can be removed if -Werror is disable. 821 for (int to{j + k}; xy != 0 && to < (2 * parts); ++to) { 822 #else 823 for (int to{j + k}; xy != 0; ++to) { 824 #endif 825 xy += product[to]; 826 product[to] = xy & partMask; 827 xy >>= partBits; 828 } 829 } 830 } 831 } 832 } 833 Integer upper{nullptr}, lower{nullptr}; 834 for (int j{0}; j < parts; ++j) { 835 lower.LEPart(j) = product[j]; 836 upper.LEPart(j) = product[j + parts]; 837 } 838 if constexpr (topPartBits < partBits) { 839 upper = upper.SHIFTL(partBits - topPartBits); 840 upper.LEPart(0) |= lower.LEPart(parts - 1) >> topPartBits; 841 lower.LEPart(parts - 1) &= topPartMask; 842 } 843 return {upper, lower}; 844 } 845 846 constexpr Product MultiplySigned(const Integer &y) const { 847 bool yIsNegative{y.IsNegative()}; 848 Integer absy{y}; 849 if (yIsNegative) { 850 absy = y.Negate().value; 851 } 852 bool isNegative{IsNegative()}; 853 Integer absx{*this}; 854 if (isNegative) { 855 absx = Negate().value; 856 } 857 Product product{absx.MultiplyUnsigned(absy)}; 858 if (isNegative != yIsNegative) { 859 product.lower = product.lower.NOT(); 860 product.upper = product.upper.NOT(); 861 Integer one{1}; 862 auto incremented{product.lower.AddUnsigned(one)}; 863 product.lower = incremented.value; 864 if (incremented.carry) { 865 product.upper = product.upper.AddUnsigned(one).value; 866 } 867 } 868 return product; 869 } 870 871 constexpr QuotientWithRemainder DivideUnsigned(const Integer &divisor) const { 872 if (divisor.IsZero()) { 873 return {MASKR(bits), Integer{}, true, false}; // overflow to max value 874 } 875 int bitsDone{LEADZ()}; 876 Integer top{SHIFTL(bitsDone)}; 877 Integer quotient, remainder; 878 for (; bitsDone < bits; ++bitsDone) { 879 auto doubledTop{top.AddUnsigned(top)}; 880 top = doubledTop.value; 881 remainder = remainder.AddUnsigned(remainder, doubledTop.carry).value; 882 bool nextBit{remainder.CompareUnsigned(divisor) != Ordering::Less}; 883 quotient = quotient.AddUnsigned(quotient, nextBit).value; 884 if (nextBit) { 885 remainder = remainder.SubtractSigned(divisor).value; 886 } 887 } 888 return {quotient, remainder, false, false}; 889 } 890 891 // A nonzero remainder has the sign of the dividend, i.e., it computes 892 // the MOD intrinsic (X-INT(X/Y)*Y), not MODULO (which is below). 893 // 8/5 = 1r3; -8/5 = -1r-3; 8/-5 = -1r3; -8/-5 = 1r-3 894 constexpr QuotientWithRemainder DivideSigned(Integer divisor) const { 895 bool dividendIsNegative{IsNegative()}; 896 bool negateQuotient{dividendIsNegative}; 897 Ordering divisorOrdering{divisor.CompareToZeroSigned()}; 898 if (divisorOrdering == Ordering::Less) { 899 negateQuotient = !negateQuotient; 900 auto negated{divisor.Negate()}; 901 if (negated.overflow) { 902 // divisor was (and is) the most negative number 903 if (CompareUnsigned(divisor) == Ordering::Equal) { 904 return {MASKR(1), Integer{}, false, bits <= 1}; 905 } else { 906 return {Integer{}, *this, false, false}; 907 } 908 } 909 divisor = negated.value; 910 } else if (divisorOrdering == Ordering::Equal) { 911 // division by zero 912 if (dividendIsNegative) { 913 return {MASKL(1), Integer{}, true, false}; 914 } else { 915 return {MASKR(bits - 1), Integer{}, true, false}; 916 } 917 } 918 Integer dividend{*this}; 919 if (dividendIsNegative) { 920 auto negated{Negate()}; 921 if (negated.overflow) { 922 // Dividend was (and remains) the most negative number. 923 // See whether the original divisor was -1 (if so, it's 1 now). 924 if (divisorOrdering == Ordering::Less && 925 divisor.CompareUnsigned(Integer{1}) == Ordering::Equal) { 926 // most negative number / -1 is the sole overflow case 927 return {*this, Integer{}, false, true}; 928 } 929 } else { 930 dividend = negated.value; 931 } 932 } 933 // Overflow is not possible, and both the dividend and divisor 934 // are now positive. 935 QuotientWithRemainder result{dividend.DivideUnsigned(divisor)}; 936 if (negateQuotient) { 937 result.quotient = result.quotient.Negate().value; 938 } 939 if (dividendIsNegative) { 940 result.remainder = result.remainder.Negate().value; 941 } 942 return result; 943 } 944 945 // Result has the sign of the divisor argument. 946 // 8 mod 5 = 3; -8 mod 5 = 2; 8 mod -5 = -2; -8 mod -5 = -3 947 constexpr ValueWithOverflow MODULO(const Integer &divisor) const { 948 bool negativeDivisor{divisor.IsNegative()}; 949 bool distinctSigns{IsNegative() != negativeDivisor}; 950 QuotientWithRemainder divided{DivideSigned(divisor)}; 951 if (distinctSigns && !divided.remainder.IsZero()) { 952 return {divided.remainder.AddUnsigned(divisor).value, divided.overflow}; 953 } else { 954 return {divided.remainder, divided.overflow}; 955 } 956 } 957 958 constexpr PowerWithErrors Power(const Integer &exponent) const { 959 PowerWithErrors result{1, false, false, false}; 960 if (exponent.IsZero()) { 961 // x**0 -> 1, including the case 0**0, which is not defined specifically 962 // in F'18 afaict; however, other Fortrans tested all produce 1, not 0, 963 // apart from nagfor, which stops with an error at runtime. 964 // Ada, APL, C's pow(), Haskell, Julia, MATLAB, and R all produce 1 too. 965 // F'77 explicitly states that 0**0 is mathematically undefined and 966 // therefore prohibited. 967 result.zeroToZero = IsZero(); 968 } else if (exponent.IsNegative()) { 969 if (IsZero()) { 970 result.divisionByZero = true; 971 result.power = MASKR(bits - 1); 972 } else if (CompareSigned(Integer{1}) == Ordering::Equal) { 973 result.power = *this; // 1**x -> 1 974 } else if (CompareSigned(Integer{-1}) == Ordering::Equal) { 975 if (exponent.BTEST(0)) { 976 result.power = *this; // (-1)**x -> -1 if x is odd 977 } 978 } else { 979 result.power.Clear(); // j**k -> 0 if |j| > 1 and k < 0 980 } 981 } else { 982 Integer shifted{*this}; 983 Integer pow{exponent}; 984 int nbits{bits - pow.LEADZ()}; 985 for (int j{0}; j < nbits; ++j) { 986 if (pow.BTEST(j)) { 987 Product product{result.power.MultiplySigned(shifted)}; 988 result.power = product.lower; 989 result.overflow |= product.SignedMultiplicationOverflowed(); 990 } 991 if (j + 1 < nbits) { 992 Product squared{shifted.MultiplySigned(shifted)}; 993 result.overflow |= squared.SignedMultiplicationOverflowed(); 994 shifted = squared.lower; 995 } 996 } 997 } 998 return result; 999 } 1000 1001 private: 1002 // A private constructor, selected by the use of nullptr, 1003 // that is used by member functions when it would be a waste 1004 // of time to initialize parts_[]. 1005 constexpr Integer(std::nullptr_t) {} 1006 1007 // Accesses parts in little-endian order. 1008 constexpr const Part &LEPart(int part) const { 1009 if constexpr (littleEndian) { 1010 return part_[part]; 1011 } else { 1012 return part_[parts - 1 - part]; 1013 } 1014 } 1015 1016 constexpr Part &LEPart(int part) { 1017 if constexpr (littleEndian) { 1018 return part_[part]; 1019 } else { 1020 return part_[parts - 1 - part]; 1021 } 1022 } 1023 1024 constexpr void SetLEPart(int part, Part x) { 1025 LEPart(part) = x & PartMask(part); 1026 } 1027 1028 static constexpr Part PartMask(int part) { 1029 return part == parts - 1 ? topPartMask : partMask; 1030 } 1031 1032 constexpr void Clear() { 1033 for (int j{0}; j < parts; ++j) { 1034 part_[j] = 0; 1035 } 1036 } 1037 1038 Part part_[parts]{}; 1039 }; 1040 1041 extern template class Integer<8>; 1042 extern template class Integer<16>; 1043 extern template class Integer<32>; 1044 extern template class Integer<64>; 1045 extern template class Integer<80>; 1046 extern template class Integer<128>; 1047 } // namespace Fortran::evaluate::value 1048 #endif // FORTRAN_EVALUATE_INTEGER_H_ 1049