564 lines · cpp
1#include "flang/Evaluate/type.h"2#include "flang/Testing/fp-testing.h"3#include "flang/Testing/testing.h"4#include "llvm/Support/raw_ostream.h"5#include <cmath>6#include <cstdio>7#include <cstdlib>8#include <type_traits>9 10using namespace Fortran::evaluate;11using namespace Fortran::common;12 13using Real2 = Scalar<Type<TypeCategory::Real, 2>>;14using Real3 = Scalar<Type<TypeCategory::Real, 3>>;15using Real4 = Scalar<Type<TypeCategory::Real, 4>>;16using Real8 = Scalar<Type<TypeCategory::Real, 8>>;17#ifdef __x86_64__18using Real10 = Scalar<Type<TypeCategory::Real, 10>>;19#endif20using Real16 = Scalar<Type<TypeCategory::Real, 16>>;21using Integer4 = Scalar<Type<TypeCategory::Integer, 4>>;22using Integer8 = Scalar<Type<TypeCategory::Integer, 8>>;23 24void dumpTest() {25 struct {26 std::uint64_t raw;27 const char *expected;28 } table[] = {29 {0x7f876543, "NaN0x7f876543"},30 {0x7f800000, "Inf"},31 {0xff800000, "-Inf"},32 {0x00000000, "0.0"},33 {0x80000000, "-0.0"},34 {0x3f800000, "0x1.0p0"},35 {0xbf800000, "-0x1.0p0"},36 {0x40000000, "0x1.0p1"},37 {0x3f000000, "0x1.0p-1"},38 {0x7f7fffff, "0x1.fffffep127"},39 {0x00800000, "0x1.0p-126"},40 {0x00400000, "0x0.8p-126"},41 {0x00000001, "0x0.000002p-126"},42 {0, nullptr},43 };44 for (int j{0}; table[j].expected != nullptr; ++j) {45 TEST(Real4{Integer4{table[j].raw}}.DumpHexadecimal() == table[j].expected)46 ("%d", j);47 }48}49 50template <typename R> void basicTests(int rm, Rounding rounding) {51 static constexpr int kind{R::bits / 8};52 char desc[64];53 using Word = typename R::Word;54 std::snprintf(desc, sizeof desc, "bits=%d, le=%d, kind=%d", R::bits,55 Word::littleEndian, kind);56 R zero;57 TEST(!zero.IsNegative())(desc);58 TEST(!zero.IsNotANumber())(desc);59 TEST(!zero.IsInfinite())(desc);60 TEST(zero.IsZero())(desc);61 MATCH(0, zero.Exponent())(desc);62 TEST(zero.RawBits().IsZero())(desc);63 MATCH(0, zero.RawBits().ToUInt64())(desc);64 TEST(zero.ABS().RawBits().IsZero())(desc);65 TEST(zero.Negate().RawBits().IEOR(Word::MASKL(1)).IsZero())(desc);66 TEST(zero.Compare(zero) == Relation::Equal)(desc);67 R minusZero{Word{std::uint64_t{1}}.SHIFTL(R::bits - 1)};68 TEST(minusZero.IsNegative())(desc);69 TEST(!minusZero.IsNotANumber())(desc);70 TEST(!minusZero.IsInfinite())(desc);71 TEST(minusZero.IsZero())(desc);72 TEST(minusZero.ABS().RawBits().IsZero())(desc);73 TEST(minusZero.Negate().RawBits().IsZero())(desc);74 MATCH(0, minusZero.Exponent())(desc);75 MATCH(0, minusZero.RawBits().LEADZ())(desc);76 MATCH(1, minusZero.RawBits().POPCNT())(desc);77 TEST(minusZero.Compare(minusZero) == Relation::Equal)(desc);78 TEST(zero.Compare(minusZero) == Relation::Equal)(desc);79 ValueWithRealFlags<R> vr;80 MATCH(0, vr.value.RawBits().ToUInt64())(desc);81 TEST(vr.flags.empty())(desc);82 R nan{Word{std::uint64_t{1}}83 .SHIFTL(R::bits)84 .SubtractSigned(Word{std::uint64_t{1}})85 .value};86 MATCH(R::bits, nan.RawBits().POPCNT())(desc);87 TEST(!nan.IsNegative())(desc);88 TEST(nan.IsNotANumber())(desc);89 TEST(!nan.IsInfinite())(desc);90 TEST(!nan.IsZero())(desc);91 TEST(zero.Compare(nan) == Relation::Unordered)(desc);92 TEST(minusZero.Compare(nan) == Relation::Unordered)(desc);93 TEST(nan.Compare(zero) == Relation::Unordered)(desc);94 TEST(nan.Compare(minusZero) == Relation::Unordered)(desc);95 TEST(nan.Compare(nan) == Relation::Unordered)(desc);96 int significandBits{R::binaryPrecision - R::isImplicitMSB};97 int exponentBits{R::bits - significandBits - 1};98 std::uint64_t maxExponent{(std::uint64_t{1} << exponentBits) - 1};99 MATCH(nan.Exponent(), maxExponent)(desc);100 Word infWord{Word{maxExponent}.SHIFTL(significandBits)};101 Word negInfWord{102 Word{maxExponent}.SHIFTL(significandBits).IOR(Word::MASKL(1))};103 if constexpr (kind == 10) { // x87104 infWord = infWord.IBSET(63);105 negInfWord = negInfWord.IBSET(63);106 }107 R inf{infWord};108 R negInf{negInfWord};109 TEST(!inf.IsNegative())(desc);110 TEST(!inf.IsNotANumber())(desc);111 TEST(inf.IsInfinite())(desc);112 TEST(!inf.IsZero())(desc);113 TEST(inf.RawBits().CompareUnsigned(inf.ABS().RawBits()) == Ordering::Equal)114 (desc);115 TEST(zero.Compare(inf) == Relation::Less)(desc);116 TEST(minusZero.Compare(inf) == Relation::Less)(desc);117 TEST(nan.Compare(inf) == Relation::Unordered)(desc);118 TEST(inf.Compare(inf) == Relation::Equal)(desc);119 TEST(negInf.IsNegative())(desc);120 TEST(!negInf.IsNotANumber())(desc);121 TEST(negInf.IsInfinite())(desc);122 TEST(!negInf.IsZero())(desc);123 TEST(inf.RawBits().CompareUnsigned(negInf.ABS().RawBits()) == Ordering::Equal)124 (desc);125 TEST(inf.RawBits().CompareUnsigned(negInf.Negate().RawBits()) ==126 Ordering::Equal)127 (desc);128 TEST(inf.Negate().RawBits().CompareUnsigned(negInf.RawBits()) ==129 Ordering::Equal)130 (desc);131 TEST(zero.Compare(negInf) == Relation::Greater)(desc);132 TEST(minusZero.Compare(negInf) == Relation::Greater)(desc);133 TEST(nan.Compare(negInf) == Relation::Unordered)(desc);134 TEST(inf.Compare(negInf) == Relation::Greater)(desc);135 TEST(negInf.Compare(negInf) == Relation::Equal)(desc);136 for (std::uint64_t j{0}; j < 63; ++j) {137 char ldesc[128];138 std::uint64_t x{1};139 x <<= j;140 std::snprintf(ldesc, sizeof ldesc, "%s j=%d x=0x%jx rm=%d", desc,141 static_cast<int>(j), static_cast<std::intmax_t>(x), rm);142 Integer8 ix{x};143 TEST(!ix.IsNegative())(ldesc);144 MATCH(x, ix.ToUInt64())(ldesc);145 vr = R::FromInteger(ix, false, rounding);146 TEST(!vr.value.IsNegative())(ldesc);147 TEST(!vr.value.IsNotANumber())(ldesc);148 TEST(!vr.value.IsZero())(ldesc);149 auto ivf = vr.value.template ToInteger<Integer8>();150 if (j > (maxExponent / 2)) {151 TEST(vr.flags.test(RealFlag::Overflow))(ldesc);152 TEST(vr.value.IsInfinite())(ldesc);153 TEST(ivf.flags.test(RealFlag::Overflow))(ldesc);154 MATCH(0x7fffffffffffffff, ivf.value.ToUInt64())(ldesc);155 } else {156 TEST(vr.flags.empty())(ldesc);157 TEST(!vr.value.IsInfinite())(ldesc);158 TEST(ivf.flags.empty())(ldesc);159 MATCH(x, ivf.value.ToUInt64())(ldesc);160 if (rounding.mode == RoundingMode::TiesToEven) { // to match stold()161 std::string decimal;162 llvm::raw_string_ostream ss{decimal};163 vr.value.AsFortran(ss, kind, false /*exact*/);164 const char *p{decimal.data()};165 MATCH(x, static_cast<std::uint64_t>(std::stold(decimal)))166 ("%s %s", ldesc, p);167 auto check{R::Read(p, rounding)};168 auto icheck{check.value.template ToInteger<Integer8>()};169 MATCH(x, icheck.value.ToUInt64())(ldesc);170 TEST(vr.value.Compare(check.value) == Relation::Equal)(ldesc);171 }172 }173 TEST(vr.value.ToWholeNumber().value.Compare(vr.value) == Relation::Equal)174 (ldesc);175 ix = ix.Negate().value;176 TEST(ix.IsNegative())(ldesc);177 x = -x;178 std::int64_t nx = x;179 MATCH(x, ix.ToUInt64())(ldesc);180 MATCH(nx, ix.ToInt64())(ldesc);181 vr = R::FromInteger(ix);182 TEST(vr.value.IsNegative())(ldesc);183 TEST(!vr.value.IsNotANumber())(ldesc);184 TEST(!vr.value.IsZero())(ldesc);185 ivf = vr.value.template ToInteger<Integer8>();186 if (j > (maxExponent / 2)) {187 TEST(vr.flags.test(RealFlag::Overflow))(ldesc);188 TEST(vr.value.IsInfinite())(ldesc);189 TEST(ivf.flags.test(RealFlag::Overflow))(ldesc);190 MATCH(0x8000000000000000, ivf.value.ToUInt64())(ldesc);191 } else {192 TEST(vr.flags.empty())(ldesc);193 TEST(!vr.value.IsInfinite())(ldesc);194 TEST(ivf.flags.empty())(ldesc);195 MATCH(x, ivf.value.ToUInt64())(ldesc);196 MATCH(nx, ivf.value.ToInt64())(ldesc);197 }198 TEST(vr.value.ToWholeNumber().value.Compare(vr.value) == Relation::Equal)199 (ldesc);200 }201}202 203// Takes an integer and distributes its bits across a floating204// point value. The LSB is used to complement the result.205std::uint32_t MakeReal(std::uint32_t n) {206 int shifts[] = {-1, 31, 23, 30, 22, 0, 24, 29, 25, 28, 26, 1, 16, 21, 2, -1};207 std::uint32_t x{0};208 for (int j{1}; shifts[j] >= 0; ++j) {209 x |= ((n >> j) & 1) << shifts[j];210 }211 x ^= -(n & 1);212 return x;213}214 215std::uint64_t MakeReal(std::uint64_t n) {216 int shifts[] = {217 -1, 63, 52, 62, 51, 0, 53, 61, 54, 60, 55, 59, 1, 16, 50, 2, -1};218 std::uint64_t x{0};219 for (int j{1}; shifts[j] >= 0; ++j) {220 x |= ((n >> j) & 1) << shifts[j];221 }222 x ^= -(n & 1);223 return x;224}225 226inline bool IsNaN(std::uint32_t x) {227 return (x & 0x7f800000) == 0x7f800000 && (x & 0x007fffff) != 0;228}229 230inline bool IsNaN(std::uint64_t x) {231 return (x & 0x7ff0000000000000) == 0x7ff0000000000000 &&232 (x & 0x000fffffffffffff) != 0;233}234 235inline bool IsInfinite(std::uint32_t x) {236 return (x & 0x7fffffff) == 0x7f800000;237}238 239inline bool IsInfinite(std::uint64_t x) {240 return (x & 0x7fffffffffffffff) == 0x7ff0000000000000;241}242 243inline bool IsNegative(std::uint32_t x) { return (x & 0x80000000) != 0; }244 245inline bool IsNegative(std::uint64_t x) {246 return (x & 0x8000000000000000) != 0;247}248 249inline std::uint32_t NormalizeNaN(std::uint32_t x) {250 if (IsNaN(x)) {251 x = 0x7fe00000;252 }253 return x;254}255 256inline std::uint64_t NormalizeNaN(std::uint64_t x) {257 if (IsNaN(x)) {258 x = 0x7ffc000000000000;259 }260 return x;261}262 263enum FlagBits {264 Overflow = 1,265 DivideByZero = 2,266 InvalidArgument = 4,267 Underflow = 8,268 Inexact = 16,269};270 271#ifdef __clang__272// clang support for fenv.h is broken, so tests of flag settings273// are disabled.274inline std::uint32_t FlagsToBits(const RealFlags &) { return 0; }275#else276inline std::uint32_t FlagsToBits(const RealFlags &flags) {277 std::uint32_t bits{0};278 if (flags.test(RealFlag::Overflow)) {279 bits |= Overflow;280 }281 if (flags.test(RealFlag::DivideByZero)) {282 bits |= DivideByZero;283 }284 if (flags.test(RealFlag::InvalidArgument)) {285 bits |= InvalidArgument;286 }287 if (flags.test(RealFlag::Underflow)) {288 bits |= Underflow;289 }290 if (flags.test(RealFlag::Inexact)) {291 bits |= Inexact;292 }293 return bits;294}295#endif // __clang__296 297template <typename UINT = std::uint32_t, typename FLT = float, typename REAL>298void inttest(std::int64_t x, int pass, Rounding rounding) {299 union {300 UINT ui;301 FLT f;302 } u;303 ScopedHostFloatingPointEnvironment fpenv;304 Integer8 ix{x};305 ValueWithRealFlags<REAL> real;306 real = real.value.FromInteger(ix, false, rounding);307#ifndef __clang__ // broken and also slow308 fpenv.ClearFlags();309#endif310 FLT fcheck = x; // TODO unsigned too311 auto actualFlags{FlagsToBits(fpenv.CurrentFlags())};312 u.f = fcheck;313 UINT rcheck{NormalizeNaN(u.ui)};314 UINT check = real.value.RawBits().ToUInt64();315 MATCH(rcheck, check)("%d 0x%llx", pass, x);316 MATCH(actualFlags, FlagsToBits(real.flags))("%d 0x%llx", pass, x);317}318 319template <typename FLT = float> FLT ToIntPower(FLT x, int power) {320 if (power == 0) {321 return x / x;322 }323 bool negative{power < 0};324 if (negative) {325 power = -power;326 }327 FLT result{1};328 while (power > 0) {329 if (power & 1) {330 result *= x;331 }332 x *= x;333 power >>= 1;334 }335 if (negative) {336 result = 1.0 / result;337 }338 return result;339}340 341template <typename FLT, int decimalDigits>342FLT TimesIntPowerOfTen(FLT x, int power) {343 if (power > decimalDigits || power < -decimalDigits) {344 auto maxExactPowerOfTen{345 TimesIntPowerOfTen<FLT, decimalDigits>(1, decimalDigits)};346 auto big{ToIntPower<FLT>(maxExactPowerOfTen, power / decimalDigits)};347 auto small{348 TimesIntPowerOfTen<FLT, decimalDigits>(1, power % decimalDigits)};349 return (x * big) * small;350 }351 return x * ToIntPower<FLT>(10.0, power);352}353 354template <typename UINT = std::uint32_t, typename FLT = float,355 typename REAL = Real4>356void subsetTests(int pass, Rounding rounding, std::uint32_t opds) {357 for (int j{0}; j < 63; ++j) {358 std::int64_t x{1};359 x <<= j;360 inttest<UINT, FLT, REAL>(x, pass, rounding);361 inttest<UINT, FLT, REAL>(-x, pass, rounding);362 }363 inttest<UINT, FLT, REAL>(0, pass, rounding);364 inttest<UINT, FLT, REAL>(365 static_cast<std::int64_t>(0x8000000000000000), pass, rounding);366 367 union {368 UINT ui;369 FLT f;370 } u;371 ScopedHostFloatingPointEnvironment fpenv;372 373 for (UINT j{0}; j < opds; ++j) {374 375 UINT rj{MakeReal(j)};376 u.ui = rj;377 FLT fj{u.f};378 REAL x{typename REAL::Word{std::uint64_t{rj}}};379 380 // unary operations381 {382 ValueWithRealFlags<REAL> aint{x.ToWholeNumber()};383#ifndef __clang__ // broken and also slow384 fpenv.ClearFlags();385#endif386 FLT fcheck{std::trunc(fj)};387 auto actualFlags{FlagsToBits(fpenv.CurrentFlags())};388 actualFlags &= ~Inexact; // x86 std::trunc can set Inexact; AINT ain't389 u.f = fcheck;390#ifndef __clang__391 if (IsNaN(u.ui)) {392 actualFlags |= InvalidArgument; // x86 std::trunc(NaN) workaround393 }394#endif395 UINT rcheck{NormalizeNaN(u.ui)};396 UINT check = aint.value.RawBits().ToUInt64();397 MATCH(rcheck, check)398 ("%d AINT(0x%jx)", pass, static_cast<std::intmax_t>(rj));399 MATCH(actualFlags, FlagsToBits(aint.flags))400 ("%d AINT(0x%jx)", pass, static_cast<std::intmax_t>(rj));401 }402 403 {404 ValueWithRealFlags<REAL> root{x.SQRT(rounding)};405#ifndef __clang__ // broken and also slow406 fpenv.ClearFlags();407#endif408 FLT fcheck{std::sqrt(fj)};409 auto actualFlags{FlagsToBits(fpenv.CurrentFlags())};410 u.f = fcheck;411 UINT rcheck{NormalizeNaN(u.ui)};412 UINT check = root.value.RawBits().ToUInt64();413 MATCH(rcheck, check)414 ("%d SQRT(0x%jx)", pass, static_cast<std::intmax_t>(rj));415 MATCH(actualFlags, FlagsToBits(root.flags))416 ("%d SQRT(0x%jx)", pass, static_cast<std::intmax_t>(rj));417 }418 419 {420 MATCH(IsNaN(rj), x.IsNotANumber())421 ("%d IsNaN(0x%jx)", pass, static_cast<std::intmax_t>(rj));422 MATCH(IsInfinite(rj), x.IsInfinite())423 ("%d IsInfinite(0x%jx)", pass, static_cast<std::intmax_t>(rj));424 425 static constexpr int kind{REAL::bits / 8};426 std::string s, cssBuf;427 llvm::raw_string_ostream ss{s};428 llvm::raw_string_ostream css{cssBuf};429 x.AsFortran(ss, kind, false /*exact*/);430 if (IsNaN(rj)) {431 css << "(0._" << kind << "/0.)";432 MATCH(cssBuf, s)433 ("%d invalid(0x%jx)", pass, static_cast<std::intmax_t>(rj));434 } else if (IsInfinite(rj)) {435 css << '(';436 if (IsNegative(rj)) {437 css << '-';438 }439 css << "1._" << kind << "/0.)";440 MATCH(cssBuf, s)441 ("%d overflow(0x%jx)", pass, static_cast<std::intmax_t>(rj));442 } else {443 const char *p = s.data();444 if (*p == '(') {445 ++p;446 }447 auto readBack{REAL::Read(p, rounding)};448 MATCH(rj, readBack.value.RawBits().ToUInt64())449 ("%d Read(AsFortran()) 0x%jx %s %g", pass,450 static_cast<std::intmax_t>(rj), s.data(), static_cast<double>(fj));451 MATCH('_', *p)452 ("%d Read(AsFortran()) 0x%jx %s %d", pass,453 static_cast<std::intmax_t>(rj), s.data(),454 static_cast<int>(p - s.data()));455 }456 }457 458 // dyadic operations459 for (UINT k{0}; k < opds; ++k) {460 UINT rk{MakeReal(k)};461 u.ui = rk;462 FLT fk{u.f};463 REAL y{typename REAL::Word{std::uint64_t{rk}}};464 {465 ValueWithRealFlags<REAL> sum{x.Add(y, rounding)};466#ifndef __clang__ // broken and also slow467 fpenv.ClearFlags();468#endif469 FLT fcheck{fj + fk};470 auto actualFlags{FlagsToBits(fpenv.CurrentFlags())};471 u.f = fcheck;472 UINT rcheck{NormalizeNaN(u.ui)};473 UINT check = sum.value.RawBits().ToUInt64();474 MATCH(rcheck, check)475 ("%d 0x%jx + 0x%jx", pass, static_cast<std::intmax_t>(rj),476 static_cast<std::intmax_t>(rk));477 MATCH(actualFlags, FlagsToBits(sum.flags))478 ("%d 0x%jx + 0x%jx", pass, static_cast<std::intmax_t>(rj),479 static_cast<std::intmax_t>(rk));480 }481 {482 ValueWithRealFlags<REAL> diff{x.Subtract(y, rounding)};483#ifndef __clang__ // broken and also slow484 fpenv.ClearFlags();485#endif486 FLT fcheck{fj - fk};487 auto actualFlags{FlagsToBits(fpenv.CurrentFlags())};488 u.f = fcheck;489 UINT rcheck{NormalizeNaN(u.ui)};490 UINT check = diff.value.RawBits().ToUInt64();491 MATCH(rcheck, check)492 ("%d 0x%jx - 0x%jx", pass, static_cast<std::intmax_t>(rj),493 static_cast<std::intmax_t>(rk));494 MATCH(actualFlags, FlagsToBits(diff.flags))495 ("%d 0x%jx - 0x%jx", pass, static_cast<std::intmax_t>(rj),496 static_cast<std::intmax_t>(rk));497 }498 {499 ValueWithRealFlags<REAL> prod{x.Multiply(y, rounding)};500#ifndef __clang__ // broken and also slow501 fpenv.ClearFlags();502#endif503 FLT fcheck{fj * fk};504 auto actualFlags{FlagsToBits(fpenv.CurrentFlags())};505 u.f = fcheck;506 UINT rcheck{NormalizeNaN(u.ui)};507 UINT check = prod.value.RawBits().ToUInt64();508 MATCH(rcheck, check)509 ("%d 0x%jx * 0x%jx", pass, static_cast<std::intmax_t>(rj),510 static_cast<std::intmax_t>(rk));511 MATCH(actualFlags, FlagsToBits(prod.flags))512 ("%d 0x%jx * 0x%jx", pass, static_cast<std::intmax_t>(rj),513 static_cast<std::intmax_t>(rk));514 }515 {516 ValueWithRealFlags<REAL> quot{x.Divide(y, rounding)};517#ifndef __clang__ // broken and also slow518 fpenv.ClearFlags();519#endif520 FLT fcheck{fj / fk};521 auto actualFlags{FlagsToBits(fpenv.CurrentFlags())};522 u.f = fcheck;523 UINT rcheck{NormalizeNaN(u.ui)};524 UINT check = quot.value.RawBits().ToUInt64();525 MATCH(rcheck, check)526 ("%d 0x%jx / 0x%jx", pass, static_cast<std::intmax_t>(rj),527 static_cast<std::intmax_t>(rk));528 MATCH(actualFlags, FlagsToBits(quot.flags))529 ("%d 0x%jx / 0x%jx", pass, static_cast<std::intmax_t>(rj),530 static_cast<std::intmax_t>(rk));531 }532 }533 }534}535 536void roundTest(int rm, Rounding rounding, std::uint32_t opds) {537 basicTests<Real2>(rm, rounding);538 basicTests<Real3>(rm, rounding);539 basicTests<Real4>(rm, rounding);540 basicTests<Real8>(rm, rounding);541#ifdef __x86_64__542 basicTests<Real10>(rm, rounding);543#endif544 basicTests<Real16>(rm, rounding);545 ScopedHostFloatingPointEnvironment::SetRounding(rounding);546 subsetTests<std::uint32_t, float, Real4>(rm, rounding, opds);547 subsetTests<std::uint64_t, double, Real8>(rm, rounding, opds);548}549 550int main() {551 dumpTest();552 std::uint32_t opds{512}; // for quick testing by default553 if (const char *p{std::getenv("REAL_TEST_OPERANDS")}) {554 // Use 8192 or 16384 for more exhaustive testing.555 opds = std::atol(p);556 }557 roundTest(0, Rounding{RoundingMode::TiesToEven}, opds);558 roundTest(1, Rounding{RoundingMode::ToZero}, opds);559 roundTest(2, Rounding{RoundingMode::Up}, opds);560 roundTest(3, Rounding{RoundingMode::Down}, opds);561 // TODO: how to test Rounding::TiesAwayFromZero on x86?562 return testing::Complete();563}564