--- /dev/null
+// Copyright (c) 2018, NVIDIA CORPORATION. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef FORTRAN_EVALUATE_FIXED_POINT_H_
+#define FORTRAN_EVALUATE_FIXED_POINT_H_
+
+// Emulates integers of a nearly arbitrary fixed size for use when the C++
+// environment does not support it. The size must be some multiple of
+// 32 bits. Signed and unsigned operations are distinct.
+
+#include "leading-zero-bit-count.h"
+#include <cinttypes>
+#include <cstddef>
+
+namespace Fortran::evaluate {
+
+enum class Ordering { Less, Equal, Greater };
+static constexpr Ordering Reverse Ordering ordering) {
+ if (ordering == Ordering::Less) {
+ return Ordering::Greater;
+ }
+ if (ordering == Ordering::Greater) {
+ return Ordering::Less;
+ }
+ return Ordering::Equal;
+}
+
+typedef <int BITS>
+class FixedPoint {
+private:
+ using Part = std::uint32_t;
+ using BigPart = std::uint64_t;
+ static constexpr int bits{BITS};
+ static constexpr int partBits{CHAR_BIT * sizeof(Part)};
+ static_assert(bits >= partBits);
+ static_assert(sizeof(BigPart) == 2 * partBits);
+ static constexpr int parts{bits / partBits};
+ static_assert(bits * partBits == parts); // no partial part
+
+public:
+ FixedPoint() = delete;
+ constexpr FixedPoint(const FixedPoint &) = default;
+ constexpr FixedPoint(std::uint64_t n) {
+ for (int j{0}; j < parts; ++j) {
+ part_[j] = n;
+ if constexpr (partBits < 64) {
+ n >>= partBits;
+ } else {
+ n = 0;
+ }
+ }
+ }
+ constexpr FixedPoint(std::int64_t n) {
+ std::int64_t signExtension{-(n < 0) << partBits};
+ for (int j{0}; j < parts; ++j) {
+ part_[j] = n;
+ if constexpr (partBits < 64) {
+ n = (n >> partBits) | signExtension;
+ } else {
+ n = signExtension;
+ }
+ }
+ }
+
+ constexpr FixedPoint &operator=(const FixedPoint &) = default;
+
+ constexpr Ordering CompareToZeroUnsigned() const {
+ for (int j{0}; j < parts; ++j) {
+ if (part_[j] != 0) {
+ return Ordering::Greater;
+ }
+ }
+ return Ordering::Equal;
+ }
+
+ constexpr Ordering CompareToZeroSigned() const {
+ if (IsNegative()) {
+ return Ordering::Less;
+ }
+ return CompareToZeroUnsigned();
+ }
+
+ constexpr Ordering CompareUnsigned(const FixedPoint &y) const {
+ for (int j{parts}; j-- > 0; ) {
+ if (part_[j] > y.part_[j]) {
+ return Ordering::Greater;
+ }
+ if (part_[j] < y.part_[j]) {
+ return Ordering::Less;
+ }
+ }
+ return Ordering::Equal;
+ }
+
+ constexpr Ordering CompareSigned(const FixedPoint &y) const {
+ if (IsNegative()) {
+ if (!y.IsNegative()) {
+ return Ordering::Less;
+ }
+ return Reverse(CompareUnsigned(y));
+ } else if (y.IsNegative()) {
+ return Ordering::Greater;
+ } else {
+ return CompareUnsigned(y);
+ }
+ }
+
+ constexpr int LeadingZeroBitCount() const {
+ for (int j{0}; j < parts; ++j) {
+ if (part_[j] != 0) {
+ return (j * partBits) + evaluate::LeadingZeroBitCount(part_[j]);
+ }
+ }
+ return bits;
+ }
+
+ constexpr std::uint64_t ToUInt64() const {
+ std::uint64_t n{0};
+ int filled{0};
+ static constexpr int toFill{bits < 64 ? bits : 64};
+ for (int j{0}; filled < 64; ++j, filled += partBits) {
+ n |= part_[j] << filled;
+ }
+ return n;
+ }
+
+ constexpr std::int64_t ToInt64() const {
+ return static_cast<std::int64_t>(ToUInt64());
+ }
+
+ constexpr void OnesComplement() {
+ for (int j{0}; j < parts; ++j) {
+ part_[j] = ~part_[j];
+ }
+ }
+
+ // Returns true on overflow (i.e., negating the most negative number)
+ constexpr bool TwosComplement() {
+ Part carry{1};
+ for (int j{0}; j < parts; ++j) {
+ Part newCarry{part_[j] == 0 && carry};
+ part_[j] = ~part_[j] + carry;
+ carry = newCarry;
+ }
+ return carry != IsNegative();
+ }
+
+ constexpr void And(const FixedPoint &y) {
+ for (int j{0}; j < parts; ++j) {
+ part_[j] &= y.part_[j];
+ }
+ }
+
+ constexpr void Or(const FixedPoint &y) {
+ for (int j{0}; j < parts; ++j) {
+ part_[j] |= y.part_[j];
+ }
+ }
+
+ constexpr void Xor(const FixedPoint &y) {
+ for (int j{0}; j < parts; ++j) {
+ part_[j] ^= y.part_[j];
+ }
+ }
+
+ constexpr void ShiftLeft(int count) {
+ if (count < 0) {
+ ShiftRight(-count);
+ } else {
+ int shiftParts{count / partBits};
+ int bitShift{count - partBits * shiftParts};
+ int j{parts-1};
+ if (bitShift == 0) {
+ for (; j >= shiftParts; --j) {
+ part_[j] = part_[j - shiftParts];
+ }
+ for (; j >= 0; --j) {
+ part_[j] = 0;
+ }
+ } else {
+ for (; j > shiftParts; --j) {
+ part_[j] = (part_[j - shiftParts] << bitShift) |
+ (part_[j - shiftParts - 1] >> (partBits - bitShift);
+ }
+ if (j == shiftParts) {
+ part_[j--] = part_[0] << bitShift;
+ }
+ for (; j >= 0; --j) {
+ part_[j] = 0;
+ }
+ }
+ }
+ }
+
+ constexpr void ShiftRightLogical(int count) { // i.e., unsigned
+ if (count < 0) {
+ ShiftLeft(-count);
+ } else {
+ int shiftParts{count / partBits};
+ int bitShift{count - partBits * shiftParts};
+ int j{0};
+ if (bitShift == 0) {
+ for (; j + shiftParts < parts; ++j) {
+ part_[j] = part_[j + shiftParts];
+ }
+ for (; j < parts; ++j) {
+ part_[j] = 0;
+ }
+ } else {
+ for (; j + shiftParts + 1 < parts; ++j) {
+ part_[j] = (part_[j + shiftParts] >> bitShift) |
+ (part_[j + shiftParts + 1] << (partBits - bitShift);
+ }
+ if (j + shiftParts + 1 == parts) {
+ part_[j++] = part_[parts - 1] >> bitShift;
+ }
+ for (; j < parts; ++j) {
+ part_[j] = 0;
+ }
+ }
+ }
+ }
+
+ // Returns carry out.
+ constexpr bool AddUnsigned(const FixedPoint &y, bool carryIn{false}) {
+ BigPart carry{carryIn};
+ for (int j{0}; j < parts; ++j) {
+ carry += part_[j];
+ part_[j] = carry += y.part_[j];
+ carry >>= 32;
+ }
+ return carry != 0;
+ }
+
+ // Returns true on overflow.
+ constexpr bool AddSigned(const FixedPoint &y) {
+ bool carry{AddUnsigned(y)};
+ return carry != IsNegative();
+ }
+
+ // Returns true on overflow.
+ constexpr bool SubtractSigned(const FixedPoint &y) {
+ FixedPoint minusy{y};
+ minusy.TwosComplement();
+ return AddSigned(minusy);
+ }
+
+ // Overwrites *this with lower half of full product.
+ constexpr void MultiplyUnsigned(const FixedPoint &y, FixedPoint &upper) {
+ Part product[2 * parts]{}; // little-endian full product
+ for (int j{0}; j < parts; ++j) {
+ if (part_[j] != 0) {
+ for (int k{0}; k < parts; ++k) {
+ if (y.part_[k] != 0) {
+ BigPart x{part_[j]};
+ x *= y.part_[k];
+ for (int to{j+k}; xy != 0; ++to) {
+ product[to] = xy += product[to];
+ xy >>= partBits;
+ }
+ }
+ }
+ }
+ }
+ for (int j{0}; j < parts; ++j) {
+ part_[j] = product[j];
+ upper.part_[j] = product[j + parts];
+ }
+ }
+
+ // Overwrites *this with lower half of full product.
+ constexpr void MultiplySigned(const FixedPoint &y, FixedPoint &upper) {
+ bool yIsNegative{y.IsNegative()};
+ FixedPoint yprime{y};
+ if (yIsNegative) {
+ yprime.TwosComplement();
+ }
+ bool isNegative{IsNegative()};
+ if (isNegative) {
+ TwosComplement();
+ }
+ MultiplyUnsigned(yprime, upper);
+ if (isNegative != yIsNegative) {
+ OnesComplement();
+ upper.OnesComplement();
+ FixedPoint one{std::uint64_t{1}};
+ if (AddUnsigned(one)) {
+ upper.AddUnsigned(one);
+ }
+ }
+ }
+
+ // Overwrites *this with quotient.
+ constexpr void DivideUnsigned(const FixedPoint &divisor, FixedPoint &remainder) {
+ FixedPoint top{*this};
+ *this = remainder = FixedPoint{0};
+ int bitsDone{top.LeadingZeroBitCount()};
+ top.ShiftLeft(bitsDone);
+ for (; bitsDone < bits; ++bitsDone) {
+ remainder.AddUnsigned(remainder, top.AddUnsigned(top));
+ bool nextBit{remainder.CompareUnsigned(divisor) != Ordering::Less};
+ quotient.AddUnsigned(quotient, nextBit);
+ if (nextBit) {
+ remainder.SubtractSigned(divisor);
+ }
+ }
+ }
+
+ // Overwrites *this with quotient. Returns true on overflow (viz.,
+ // the most negative value divided by -1) or division by zero.
+ constexpr bool DivideSigned(FixedPoint divisor, FixedPoint &remainder) {
+ bool negateQuotient{false}, negateRemainder{false};
+ if (IsNegative()) {
+ negateQuotient = negateRemainder = true;
+ TwosComplement();
+ }
+ Ordering divisorOrdering{divisor.CompareToZeroSigned()};
+ bool overflow{divisorOrdering == Ordering::Equal};
+ if (divisorOrdering == Ordering::Less) {
+ negateQuotient = !negateQuotient;
+ divisor.TwosComplement();
+ }
+ DivideUnsigned(divisor, remainder);
+ overflow |= IsNegative();
+ if (negateQuotient) {
+ TwosComplement();
+ }
+ if (negateRemainder) {
+ remainder.TwosComplement();
+ }
+ return overflow;
+ }
+
+private:
+ constexpr bool IsNegative() const {
+ return (part_[parts-1] >> (partBits - 1)) & 1;
+ }
+
+ Part part_[parts]; // little-endian order: [parts-1] is most significant
+};
+} // namespace Fortran::evaluate
+#endif // FORTRAN_EVALUATE_FIXED_POINT_H_
--- /dev/null
+// Copyright (c) 2018, NVIDIA CORPORATION. All rights reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef FORTRAN_EVALUATE_LEADING_ZERO_BIT_COUNT_H_
+#define FORTRAN_EVALUATE_LEADING_ZERO_BIT_COUNT_H_
+
+// A fast and portable function that counts the number of leading zero bits
+// in an integer value. (If the most significant bit is set, the leading
+// zero count is zero; if no bit is set, the leading zero count is the
+// word size in bits; otherwise, it's the largest left shift count that
+// doesn't reduce the number of bits in the word that are set.)
+
+#include <cinttypes>
+
+namespace Fortran::evaluate {
+namespace {
+// The following magic constant is a binary deBruijn sequence.
+// It has the remarkable property that if one extends it
+// (virtually) on the right with 5 more zero bits, then all
+// of the 64 contiguous framed blocks of six bits in the
+// extended 69-bit sequence are distinct. Consequently,
+// if one shifts it left by any shift count [0..63] with
+// truncation and extracts the uppermost six bit field
+// of the shifted value, each shift count maps to a distinct
+// field value. That means that we can map those 64 field
+// values back to the shift counts that produce them,
+// and (the point) this means that we can shift this value
+// by an unknown bit count in [0..63] and then figure out
+// what that count must have been.
+// 0 7 e d d 5 e 5 9 a 4 e 2 8 c 2
+// 0000011111101101110101011110010110011010010011100010100011000010
+static constexpr std::uint64_t deBruijn{0x07edd5e59a4e28c2};
+static constexpr std::uint8_t mapping[64]{
+ 63, 0, 58, 1, 59, 47, 53, 2, 60, 39, 48, 27, 54, 33, 42, 3,
+ 61, 51, 37, 40, 49, 18, 28, 20, 55, 30, 34, 11, 43, 14, 22, 4,
+ 62, 57, 46, 52, 38, 26, 32, 41, 50, 36, 17, 19, 29, 10, 13, 21,
+ 56, 45, 25, 31, 35, 16, 9, 12, 44, 24, 15, 8, 23, 7, 6, 5 };
+} // namespace
+
+inline constexpr int LeadingZeroBitCount(std::uint64_t x) {
+ if (x == 0) {
+ return 64;
+ } else {
+ x |= x >> 1;
+ x |= x >> 2;
+ x |= x >> 4;
+ x |= x >> 8;
+ x |= x >> 16;
+ x |= x >> 32;
+ // All of the bits below the uppermost set bit are now also set.
+ x -= x >> 1; // All of the bits below the uppermost are now clear.
+ // x now has exactly one bit set, so it is a power of two, so
+ // multiplication by x is equivalent to a left shift by its
+ // base-2 logarithm. We calculate that unknown base-2 logarithm
+ // by shifting the deBruijn sequence and mapping the framed value.
+ int base2Log{mapping[(x * deBruijn) >> 58]};
+ return 63 - base2Log; // convert to leading zero count
+ }
+}
+
+inline constexpr int LeadingZeroBitCount(std::uint32_t x) {
+ return LeadingZeroBitCount(static_cast<std::uint64_t>(x)) - 32;
+}
+
+inline constexpr int LeadingZeroBitCount(std::uint16_t x) {
+ return LeadingZeroBitCount(static_cast<std::uint64_t>(x)) - 48;
+}
+
+namespace {
+static constexpr std::uint8_t eightBitLeadingZeroBitCount[256]{
+ 8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
+ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
+ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
+ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
+ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+};
+} // namespace
+
+inline constexpr int LeadingZeroBitCount(std::uint8_t x) {
+ return eightBitLeadingZeroBitCount[x];
+}
+} // namespace Fortran::evaluate
+#endif // FORTRAN_EVALUATE_LEADING_ZERO_BIT_COUNT_H_