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 CompareUnsigned(that) == Ordering::Less; 180 } 181 constexpr bool operator<=(const Integer &that) const { 182 return CompareUnsigned(that) != Ordering::Greater; 183 } 184 constexpr bool operator==(const Integer &that) const { 185 return CompareUnsigned(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 CompareUnsigned(that) != Ordering::Less; 192 } 193 constexpr bool operator>(const Integer &that) const { 194 return CompareUnsigned(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); } 361 static constexpr int RANGE{// in the sense of SELECTED_INT_KIND 362 // This magic value is LOG10(2.)*1E12. 363 static_cast<int>(((bits - 1) * 301029995664) / 1000000000000)}; 364 IsZero()365 constexpr bool IsZero() const { 366 for (int j{0}; j < parts; ++j) { 367 if (part_[j] != 0) { 368 return false; 369 } 370 } 371 return true; 372 } 373 IsNegative()374 constexpr bool IsNegative() const { 375 return (LEPart(parts - 1) >> (topPartBits - 1)) & 1; 376 } 377 CompareToZeroSigned()378 constexpr Ordering CompareToZeroSigned() const { 379 if (IsNegative()) { 380 return Ordering::Less; 381 } else if (IsZero()) { 382 return Ordering::Equal; 383 } else { 384 return Ordering::Greater; 385 } 386 } 387 388 // Count the number of contiguous most-significant bit positions 389 // that are clear. LEADZ()390 constexpr int LEADZ() const { 391 if (LEPart(parts - 1) != 0) { 392 int lzbc{common::LeadingZeroBitCount(LEPart(parts - 1))}; 393 return lzbc - extraTopPartBits; 394 } 395 int upperZeroes{topPartBits}; 396 for (int j{1}; j < parts; ++j) { 397 if (Part p{LEPart(parts - 1 - j)}) { 398 int lzbc{common::LeadingZeroBitCount(p)}; 399 return upperZeroes + lzbc - extraPartBits; 400 } 401 upperZeroes += partBits; 402 } 403 return bits; 404 } 405 406 // Count the number of bit positions that are set. POPCNT()407 constexpr int POPCNT() const { 408 int count{0}; 409 for (int j{0}; j < parts; ++j) { 410 count += common::BitPopulationCount(part_[j]); 411 } 412 return count; 413 } 414 415 // True when POPCNT is odd. POPPAR()416 constexpr bool POPPAR() const { return POPCNT() & 1; } 417 TRAILZ()418 constexpr int TRAILZ() const { 419 auto minus1{AddUnsigned(MASKR(bits))}; // { x-1, carry = x > 0 } 420 if (!minus1.carry) { 421 return bits; // was zero 422 } else { 423 // x ^ (x-1) has all bits set at and below original least-order set bit. 424 return IEOR(minus1.value).POPCNT() - 1; 425 } 426 } 427 BTEST(int pos)428 constexpr bool BTEST(int pos) const { 429 if (pos < 0 || pos >= bits) { 430 return false; 431 } else { 432 return (LEPart(pos / partBits) >> (pos % partBits)) & 1; 433 } 434 } 435 CompareUnsigned(const Integer & y)436 constexpr Ordering CompareUnsigned(const Integer &y) const { 437 for (int j{parts}; j-- > 0;) { 438 if (LEPart(j) > y.LEPart(j)) { 439 return Ordering::Greater; 440 } 441 if (LEPart(j) < y.LEPart(j)) { 442 return Ordering::Less; 443 } 444 } 445 return Ordering::Equal; 446 } 447 BGE(const Integer & y)448 constexpr bool BGE(const Integer &y) const { 449 return CompareUnsigned(y) != Ordering::Less; 450 } BGT(const Integer & y)451 constexpr bool BGT(const Integer &y) const { 452 return CompareUnsigned(y) == Ordering::Greater; 453 } BLE(const Integer & y)454 constexpr bool BLE(const Integer &y) const { return !BGT(y); } BLT(const Integer & y)455 constexpr bool BLT(const Integer &y) const { return !BGE(y); } 456 CompareSigned(const Integer & y)457 constexpr Ordering CompareSigned(const Integer &y) const { 458 bool isNegative{IsNegative()}; 459 if (isNegative != y.IsNegative()) { 460 return isNegative ? Ordering::Less : Ordering::Greater; 461 } 462 return CompareUnsigned(y); 463 } 464 ToUInt64()465 constexpr std::uint64_t ToUInt64() const { 466 std::uint64_t n{LEPart(0)}; 467 int filled{partBits}; 468 for (int j{1}; filled < 64 && j < parts; ++j, filled += partBits) { 469 n |= std::uint64_t{LEPart(j)} << filled; 470 } 471 return n; 472 } 473 ToInt64()474 constexpr std::int64_t ToInt64() const { 475 std::int64_t signExtended = ToUInt64(); 476 if constexpr (bits < 64) { 477 signExtended |= -(signExtended >> (bits - 1)) << bits; 478 } 479 return signExtended; 480 } 481 482 // Ones'-complement (i.e., C's ~) NOT()483 constexpr Integer NOT() const { 484 Integer result{nullptr}; 485 for (int j{0}; j < parts; ++j) { 486 result.SetLEPart(j, ~LEPart(j)); 487 } 488 return result; 489 } 490 491 // Two's-complement negation (-x = ~x + 1). 492 // An overflow flag accompanies the result, and will be true when the 493 // operand is the most negative signed number (MASKL(1)). Negate()494 constexpr ValueWithOverflow Negate() const { 495 Integer result{nullptr}; 496 Part carry{1}; 497 for (int j{0}; j + 1 < parts; ++j) { 498 Part newCarry{LEPart(j) == 0 && carry}; 499 result.SetLEPart(j, ~LEPart(j) + carry); 500 carry = newCarry; 501 } 502 Part top{LEPart(parts - 1)}; 503 result.SetLEPart(parts - 1, ~top + carry); 504 bool overflow{top != 0 && result.LEPart(parts - 1) == top}; 505 return {result, overflow}; 506 } 507 ABS()508 constexpr ValueWithOverflow ABS() const { 509 if (IsNegative()) { 510 return Negate(); 511 } else { 512 return {*this, false}; 513 } 514 } 515 516 // Shifts the operand left when the count is positive, right when negative. 517 // Vacated bit positions are filled with zeroes. ISHFT(int count)518 constexpr Integer ISHFT(int count) const { 519 if (count < 0) { 520 return SHIFTR(-count); 521 } else { 522 return SHIFTL(count); 523 } 524 } 525 526 // Left shift with zero fill. SHIFTL(int count)527 constexpr Integer SHIFTL(int count) const { 528 if (count <= 0) { 529 return *this; 530 } else { 531 Integer result{nullptr}; 532 int shiftParts{count / partBits}; 533 int bitShift{count - partBits * shiftParts}; 534 int j{parts - 1}; 535 if (bitShift == 0) { 536 for (; j >= shiftParts; --j) { 537 result.SetLEPart(j, LEPart(j - shiftParts)); 538 } 539 for (; j >= 0; --j) { 540 result.LEPart(j) = 0; 541 } 542 } else { 543 for (; j > shiftParts; --j) { 544 result.SetLEPart(j, 545 ((LEPart(j - shiftParts) << bitShift) | 546 (LEPart(j - shiftParts - 1) >> (partBits - bitShift)))); 547 } 548 if (j == shiftParts) { 549 result.SetLEPart(j, LEPart(0) << bitShift); 550 --j; 551 } 552 for (; j >= 0; --j) { 553 result.LEPart(j) = 0; 554 } 555 } 556 return result; 557 } 558 } 559 560 // Circular shift of a field of least-significant bits. The least-order 561 // "size" bits are shifted circularly in place by "count" positions; 562 // the shift is leftward if count is nonnegative, rightward otherwise. 563 // Higher-order bits are unchanged. 564 constexpr Integer ISHFTC(int count, int size = bits) const { 565 if (count == 0 || size <= 0) { 566 return *this; 567 } 568 if (size > bits) { 569 size = bits; 570 } 571 count %= size; 572 if (count == 0) { 573 return *this; 574 } 575 int middleBits{size - count}, leastBits{count}; 576 if (count < 0) { 577 middleBits = -count; 578 leastBits = size + count; 579 } 580 if (size == bits) { 581 return SHIFTL(leastBits).IOR(SHIFTR(middleBits)); 582 } 583 Integer unchanged{IAND(MASKL(bits - size))}; 584 Integer middle{IAND(MASKR(middleBits)).SHIFTL(leastBits)}; 585 Integer least{SHIFTR(middleBits).IAND(MASKR(leastBits))}; 586 return unchanged.IOR(middle).IOR(least); 587 } 588 589 // Double shifts, aka shifts with specific fill. SHIFTLWithFill(const Integer & fill,int count)590 constexpr Integer SHIFTLWithFill(const Integer &fill, int count) const { 591 if (count <= 0) { 592 return *this; 593 } else if (count >= 2 * bits) { 594 return {}; 595 } else if (count > bits) { 596 return fill.SHIFTL(count - bits); 597 } else if (count == bits) { 598 return fill; 599 } else { 600 return SHIFTL(count).IOR(fill.SHIFTR(bits - count)); 601 } 602 } 603 SHIFTRWithFill(const Integer & fill,int count)604 constexpr Integer SHIFTRWithFill(const Integer &fill, int count) const { 605 if (count <= 0) { 606 return *this; 607 } else if (count >= 2 * bits) { 608 return {}; 609 } else if (count > bits) { 610 return fill.SHIFTR(count - bits); 611 } else if (count == bits) { 612 return fill; 613 } else { 614 return SHIFTR(count).IOR(fill.SHIFTL(bits - count)); 615 } 616 } 617 DSHIFTL(const Integer & fill,int count)618 constexpr Integer DSHIFTL(const Integer &fill, int count) const { 619 // DSHIFTL(I,J) shifts I:J left; the second argument is the right fill. 620 return SHIFTLWithFill(fill, count); 621 } 622 DSHIFTR(const Integer & value,int count)623 constexpr Integer DSHIFTR(const Integer &value, int count) const { 624 // DSHIFTR(I,J) shifts I:J right; the *first* argument is the left fill. 625 return value.SHIFTRWithFill(*this, count); 626 } 627 628 // Vacated upper bits are filled with zeroes. SHIFTR(int count)629 constexpr Integer SHIFTR(int count) const { 630 if (count <= 0) { 631 return *this; 632 } else { 633 Integer result{nullptr}; 634 int shiftParts{count / partBits}; 635 int bitShift{count - partBits * shiftParts}; 636 int j{0}; 637 if (bitShift == 0) { 638 for (; j + shiftParts < parts; ++j) { 639 result.LEPart(j) = LEPart(j + shiftParts); 640 } 641 for (; j < parts; ++j) { 642 result.LEPart(j) = 0; 643 } 644 } else { 645 for (; j + shiftParts + 1 < parts; ++j) { 646 result.SetLEPart(j, 647 (LEPart(j + shiftParts) >> bitShift) | 648 (LEPart(j + shiftParts + 1) << (partBits - bitShift))); 649 } 650 if (j + shiftParts + 1 == parts) { 651 result.LEPart(j++) = LEPart(parts - 1) >> bitShift; 652 } 653 for (; j < parts; ++j) { 654 result.LEPart(j) = 0; 655 } 656 } 657 return result; 658 } 659 } 660 661 // Be advised, an arithmetic (sign-filling) right shift is not 662 // the same as a division by a power of two in all cases. SHIFTA(int count)663 constexpr Integer SHIFTA(int count) const { 664 if (count <= 0) { 665 return *this; 666 } else if (IsNegative()) { 667 return SHIFTR(count).IOR(MASKL(count)); 668 } else { 669 return SHIFTR(count); 670 } 671 } 672 673 // Clears a single bit. IBCLR(int pos)674 constexpr Integer IBCLR(int pos) const { 675 if (pos < 0 || pos >= bits) { 676 return *this; 677 } else { 678 Integer result{*this}; 679 result.LEPart(pos / partBits) &= ~(Part{1} << (pos % partBits)); 680 return result; 681 } 682 } 683 684 // Sets a single bit. IBSET(int pos)685 constexpr Integer IBSET(int pos) const { 686 if (pos < 0 || pos >= bits) { 687 return *this; 688 } else { 689 Integer result{*this}; 690 result.LEPart(pos / partBits) |= Part{1} << (pos % partBits); 691 return result; 692 } 693 } 694 695 // Extracts a field. IBITS(int pos,int size)696 constexpr Integer IBITS(int pos, int size) const { 697 return SHIFTR(pos).IAND(MASKR(size)); 698 } 699 IAND(const Integer & y)700 constexpr Integer IAND(const Integer &y) const { 701 Integer result{nullptr}; 702 for (int j{0}; j < parts; ++j) { 703 result.LEPart(j) = LEPart(j) & y.LEPart(j); 704 } 705 return result; 706 } 707 IOR(const Integer & y)708 constexpr Integer IOR(const Integer &y) const { 709 Integer result{nullptr}; 710 for (int j{0}; j < parts; ++j) { 711 result.LEPart(j) = LEPart(j) | y.LEPart(j); 712 } 713 return result; 714 } 715 IEOR(const Integer & y)716 constexpr Integer IEOR(const Integer &y) const { 717 Integer result{nullptr}; 718 for (int j{0}; j < parts; ++j) { 719 result.LEPart(j) = LEPart(j) ^ y.LEPart(j); 720 } 721 return result; 722 } 723 MERGE_BITS(const Integer & y,const Integer & mask)724 constexpr Integer MERGE_BITS(const Integer &y, const Integer &mask) const { 725 return IAND(mask).IOR(y.IAND(mask.NOT())); 726 } 727 MAX(const Integer & y)728 constexpr Integer MAX(const Integer &y) const { 729 if (CompareSigned(y) == Ordering::Less) { 730 return y; 731 } else { 732 return *this; 733 } 734 } 735 MIN(const Integer & y)736 constexpr Integer MIN(const Integer &y) const { 737 if (CompareSigned(y) == Ordering::Less) { 738 return *this; 739 } else { 740 return y; 741 } 742 } 743 744 // Unsigned addition with carry. 745 constexpr ValueWithCarry AddUnsigned( 746 const Integer &y, bool carryIn = false) const { 747 Integer sum{nullptr}; 748 BigPart carry{carryIn}; 749 for (int j{0}; j + 1 < parts; ++j) { 750 carry += LEPart(j); 751 carry += y.LEPart(j); 752 sum.SetLEPart(j, carry); 753 carry >>= partBits; 754 } 755 carry += LEPart(parts - 1); 756 carry += y.LEPart(parts - 1); 757 sum.SetLEPart(parts - 1, carry); 758 return {sum, carry > topPartMask}; 759 } 760 AddSigned(const Integer & y)761 constexpr ValueWithOverflow AddSigned(const Integer &y) const { 762 bool isNegative{IsNegative()}; 763 bool sameSign{isNegative == y.IsNegative()}; 764 ValueWithCarry sum{AddUnsigned(y)}; 765 bool overflow{sameSign && sum.value.IsNegative() != isNegative}; 766 return {sum.value, overflow}; 767 } 768 SubtractSigned(const Integer & y)769 constexpr ValueWithOverflow SubtractSigned(const Integer &y) const { 770 bool isNegative{IsNegative()}; 771 bool sameSign{isNegative == y.IsNegative()}; 772 ValueWithCarry diff{AddUnsigned(y.Negate().value)}; 773 bool overflow{!sameSign && diff.value.IsNegative() != isNegative}; 774 return {diff.value, overflow}; 775 } 776 777 // MAX(X-Y, 0) DIM(const Integer & y)778 constexpr Integer DIM(const Integer &y) const { 779 if (CompareSigned(y) != Ordering::Greater) { 780 return {}; 781 } else { 782 return SubtractSigned(y).value; 783 } 784 } 785 SIGN(bool toNegative)786 constexpr ValueWithOverflow SIGN(bool toNegative) const { 787 if (toNegative == IsNegative()) { 788 return {*this, false}; 789 } else if (toNegative) { 790 return Negate(); 791 } else { 792 return ABS(); 793 } 794 } 795 SIGN(const Integer & sign)796 constexpr ValueWithOverflow SIGN(const Integer &sign) const { 797 return SIGN(sign.IsNegative()); 798 } 799 MultiplyUnsigned(const Integer & y)800 constexpr Product MultiplyUnsigned(const Integer &y) const { 801 Part product[2 * parts]{}; // little-endian full product 802 for (int j{0}; j < parts; ++j) { 803 if (Part xpart{LEPart(j)}) { 804 for (int k{0}; k < parts; ++k) { 805 if (Part ypart{y.LEPart(k)}) { 806 BigPart xy{xpart}; 807 xy *= ypart; 808 #if defined __GNUC__ && __GNUC__ < 8 809 // && to < (2 * parts) was added to avoid GCC < 8 build failure on 810 // -Werror=array-bounds. This can be removed if -Werror is disable. 811 for (int to{j + k}; xy != 0 && to < (2 * parts); ++to) { 812 #else 813 for (int to{j + k}; xy != 0; ++to) { 814 #endif 815 xy += product[to]; 816 product[to] = xy & partMask; 817 xy >>= partBits; 818 } 819 } 820 } 821 } 822 } 823 Integer upper{nullptr}, lower{nullptr}; 824 for (int j{0}; j < parts; ++j) { 825 lower.LEPart(j) = product[j]; 826 upper.LEPart(j) = product[j + parts]; 827 } 828 if constexpr (topPartBits < partBits) { 829 upper = upper.SHIFTL(partBits - topPartBits); 830 upper.LEPart(0) |= lower.LEPart(parts - 1) >> topPartBits; 831 lower.LEPart(parts - 1) &= topPartMask; 832 } 833 return {upper, lower}; 834 } 835 836 constexpr Product MultiplySigned(const Integer &y) const { 837 bool yIsNegative{y.IsNegative()}; 838 Integer absy{y}; 839 if (yIsNegative) { 840 absy = y.Negate().value; 841 } 842 bool isNegative{IsNegative()}; 843 Integer absx{*this}; 844 if (isNegative) { 845 absx = Negate().value; 846 } 847 Product product{absx.MultiplyUnsigned(absy)}; 848 if (isNegative != yIsNegative) { 849 product.lower = product.lower.NOT(); 850 product.upper = product.upper.NOT(); 851 Integer one{1}; 852 auto incremented{product.lower.AddUnsigned(one)}; 853 product.lower = incremented.value; 854 if (incremented.carry) { 855 product.upper = product.upper.AddUnsigned(one).value; 856 } 857 } 858 return product; 859 } 860 861 constexpr QuotientWithRemainder DivideUnsigned(const Integer &divisor) const { 862 if (divisor.IsZero()) { 863 return {MASKR(bits), Integer{}, true, false}; // overflow to max value 864 } 865 int bitsDone{LEADZ()}; 866 Integer top{SHIFTL(bitsDone)}; 867 Integer quotient, remainder; 868 for (; bitsDone < bits; ++bitsDone) { 869 auto doubledTop{top.AddUnsigned(top)}; 870 top = doubledTop.value; 871 remainder = remainder.AddUnsigned(remainder, doubledTop.carry).value; 872 bool nextBit{remainder.CompareUnsigned(divisor) != Ordering::Less}; 873 quotient = quotient.AddUnsigned(quotient, nextBit).value; 874 if (nextBit) { 875 remainder = remainder.SubtractSigned(divisor).value; 876 } 877 } 878 return {quotient, remainder, false, false}; 879 } 880 881 // A nonzero remainder has the sign of the dividend, i.e., it computes 882 // the MOD intrinsic (X-INT(X/Y)*Y), not MODULO (which is below). 883 // 8/5 = 1r3; -8/5 = -1r-3; 8/-5 = -1r3; -8/-5 = 1r-3 884 constexpr QuotientWithRemainder DivideSigned(Integer divisor) const { 885 bool dividendIsNegative{IsNegative()}; 886 bool negateQuotient{dividendIsNegative}; 887 Ordering divisorOrdering{divisor.CompareToZeroSigned()}; 888 if (divisorOrdering == Ordering::Less) { 889 negateQuotient = !negateQuotient; 890 auto negated{divisor.Negate()}; 891 if (negated.overflow) { 892 // divisor was (and is) the most negative number 893 if (CompareUnsigned(divisor) == Ordering::Equal) { 894 return {MASKR(1), Integer{}, false, bits <= 1}; 895 } else { 896 return {Integer{}, *this, false, false}; 897 } 898 } 899 divisor = negated.value; 900 } else if (divisorOrdering == Ordering::Equal) { 901 // division by zero 902 if (dividendIsNegative) { 903 return {MASKL(1), Integer{}, true, false}; 904 } else { 905 return {MASKR(bits - 1), Integer{}, true, false}; 906 } 907 } 908 Integer dividend{*this}; 909 if (dividendIsNegative) { 910 auto negated{Negate()}; 911 if (negated.overflow) { 912 // Dividend was (and remains) the most negative number. 913 // See whether the original divisor was -1 (if so, it's 1 now). 914 if (divisorOrdering == Ordering::Less && 915 divisor.CompareUnsigned(Integer{1}) == Ordering::Equal) { 916 // most negative number / -1 is the sole overflow case 917 return {*this, Integer{}, false, true}; 918 } 919 } else { 920 dividend = negated.value; 921 } 922 } 923 // Overflow is not possible, and both the dividend and divisor 924 // are now positive. 925 QuotientWithRemainder result{dividend.DivideUnsigned(divisor)}; 926 if (negateQuotient) { 927 result.quotient = result.quotient.Negate().value; 928 } 929 if (dividendIsNegative) { 930 result.remainder = result.remainder.Negate().value; 931 } 932 return result; 933 } 934 935 // Result has the sign of the divisor argument. 936 // 8 mod 5 = 3; -8 mod 5 = 2; 8 mod -5 = -2; -8 mod -5 = -3 937 constexpr ValueWithOverflow MODULO(const Integer &divisor) const { 938 bool negativeDivisor{divisor.IsNegative()}; 939 bool distinctSigns{IsNegative() != negativeDivisor}; 940 QuotientWithRemainder divided{DivideSigned(divisor)}; 941 if (distinctSigns && !divided.remainder.IsZero()) { 942 return {divided.remainder.AddUnsigned(divisor).value, divided.overflow}; 943 } else { 944 return {divided.remainder, divided.overflow}; 945 } 946 } 947 948 constexpr PowerWithErrors Power(const Integer &exponent) const { 949 PowerWithErrors result{1, false, false, false}; 950 if (exponent.IsZero()) { 951 // x**0 -> 1, including the case 0**0, which is not defined specifically 952 // in F'18 afaict; however, other Fortrans tested all produce 1, not 0, 953 // apart from nagfor, which stops with an error at runtime. 954 // Ada, APL, C's pow(), Haskell, Julia, MATLAB, and R all produce 1 too. 955 // F'77 explicitly states that 0**0 is mathematically undefined and 956 // therefore prohibited. 957 result.zeroToZero = IsZero(); 958 } else if (exponent.IsNegative()) { 959 if (IsZero()) { 960 result.divisionByZero = true; 961 result.power = MASKR(bits - 1); 962 } else if (CompareSigned(Integer{1}) == Ordering::Equal) { 963 result.power = *this; // 1**x -> 1 964 } else if (CompareSigned(Integer{-1}) == Ordering::Equal) { 965 if (exponent.BTEST(0)) { 966 result.power = *this; // (-1)**x -> -1 if x is odd 967 } 968 } else { 969 result.power.Clear(); // j**k -> 0 if |j| > 1 and k < 0 970 } 971 } else { 972 Integer shifted{*this}; 973 Integer pow{exponent}; 974 int nbits{bits - pow.LEADZ()}; 975 for (int j{0}; j < nbits; ++j) { 976 if (pow.BTEST(j)) { 977 Product product{result.power.MultiplySigned(shifted)}; 978 result.power = product.lower; 979 result.overflow |= product.SignedMultiplicationOverflowed(); 980 } 981 if (j + 1 < nbits) { 982 Product squared{shifted.MultiplySigned(shifted)}; 983 result.overflow |= squared.SignedMultiplicationOverflowed(); 984 shifted = squared.lower; 985 } 986 } 987 } 988 return result; 989 } 990 991 private: 992 // A private constructor, selected by the use of nullptr, 993 // that is used by member functions when it would be a waste 994 // of time to initialize parts_[]. 995 constexpr Integer(std::nullptr_t) {} 996 997 // Accesses parts in little-endian order. 998 constexpr const Part &LEPart(int part) const { 999 if constexpr (littleEndian) { 1000 return part_[part]; 1001 } else { 1002 return part_[parts - 1 - part]; 1003 } 1004 } 1005 1006 constexpr Part &LEPart(int part) { 1007 if constexpr (littleEndian) { 1008 return part_[part]; 1009 } else { 1010 return part_[parts - 1 - part]; 1011 } 1012 } 1013 1014 constexpr void SetLEPart(int part, Part x) { 1015 LEPart(part) = x & PartMask(part); 1016 } 1017 1018 static constexpr Part PartMask(int part) { 1019 return part == parts - 1 ? topPartMask : partMask; 1020 } 1021 1022 constexpr void Clear() { 1023 for (int j{0}; j < parts; ++j) { 1024 part_[j] = 0; 1025 } 1026 } 1027 1028 Part part_[parts]{}; 1029 }; 1030 1031 extern template class Integer<8>; 1032 extern template class Integer<16>; 1033 extern template class Integer<32>; 1034 extern template class Integer<64>; 1035 extern template class Integer<80>; 1036 extern template class Integer<128>; 1037 } // namespace Fortran::evaluate::value 1038 #endif // FORTRAN_EVALUATE_INTEGER_H_ 1039